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Abstract 

The article reviews the current status of a theoretical approach to the problem of the 
emission of gravitational waves by isolated systems in the context of general relativity. Part [I] 
of the article deals with general post-Newtonian sources. The exterior field of the source is 
investigated by means of a combination of analytic post-Minkowskian and multipolar approx- 
imations. The physical observables in the far-zone of the source are described by a specific set 
of radiative multipole moments. By matching the exterior solution to the metric of the post- 
Newtonian source in the near-zone we obtain the explicit expressions of the source multipole 
moments. The relationships between the radiative and source moments involve many non- 
linear multipole interactions, among them those associated with the tails (and tails-of-tails) 
of gravitational waves. Part [n] of the article is devoted to the application to compact bi- 
nary systems. We present the equations of binary motion, and the associated Lagrangian and 
Hamiltonian, at the third post-Newtonian (3PN) order beyond the Newtonian acceleration. 
The gravitational-wave energy flux, taking consistently into account the relativistic corrections 
in the binary moments as well as the various tail effects, is derived through 3.5PN order with 
respect to the quadrupole formalism. The binary's orbital phase, whose prior knowledge is 
crucial for searching and analyzing the signals from inspiralling compact binaries, is deduced 
from an energy balance argument. 
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1 Introduction 



The theory of gravitational radiation from isolated sources, in the context of general relativity, is 
a fascinating science that can be explored by means of what was referred to in the French XVIIIth 
century as I 'analyse sublime: the analytical (i.e. mathematical) method, and more specifically 
the resolution of partial differential equations. Indeed, the field equations of general relativity, 
when use is made of the harmonic-coordinate conditions, take the form of a quasi-linear hyperbolic 
differential system of equations, involving the famous wave operator or d'Alembertian (denoted 
□), invented by d'Alembert in his Traite de dynamique of 1743. 

Nowadays, the importance of the field lies in the exciting possibility of comparing the theory 
with contemporary astrophysical observations, made by a new generation of detectors - large- 
scale optical interferometers LIGO, VIRGO, GEO and TAMA - that should routinely observe the 
gravitational waves produced by massive and rapidly evolving systems such as inspiralling compact 
binaries. To prepare these experiments, the required theoretical work consists of carrying out a 
sufficiently general solution of the Einstein field equations, valid for a large class of matter systems, 
and describing the physical processes of the emission and propagation of the waves from the source 
to the distant detector, as well as their back-reaction onto the source. 

1.1 Gravitational- wave generation formalisms 

The basic problem we face is to relate the asymptotic gravitational-wave form hij generated by 
some isolated source, at the location of some detector in the wave zone of the source, to the stress- 
energy tensor T a @ of the matter fields 1 . For general sources it is hopeless to solve the problem 
via a rigorous deduction within the exact theory of general relativity, and we have to resort to 
approximation methods, keeping in mind that, sadly, such methods are often not related in a very 
precise mathematical way to the first principles of the theory. Therefore, a general wave-generation 
formalism must solve the field equations, and the non-linearity therein, by imposing some suitable 
approximation series in one or several small physical parameters. Of ourse the ultimate aim of 
approximation methods is to extract from the theory some firm predictions for the outcome of 
experiments such as VIRGO and LIGO. Some important approximations that we shall use in 
this article are the post-Newtonian method (or non-linear 1/c-expansion), the post-Minkowskian 
method or non-linear iteration (G-expansion), the multipole decomposition in irreducible represen- 
tations of the rotation group (or equivalently et-expansion in the source radius) , and the far-zone 
expansion (l/i?-cxpansion in the distance). In particular, the post-Newtonian expansion has pro- 
vided us in the past with our best insights into the problems of motion and radiation in general 
relativity. The most successful wave-generation formalisms make a gourmet cocktail of all these 
approximation methods. For reviews on analytic approximations and applications to the motion 
and the gravitational wave-generation see Refs. [271%! l%3l IM1 127131 12151 IP71 122] . 

The post-Newtonian approximation is valid under the assumptions of a weak gravitational field 
inside the source (we shall see later how to model neutron stars and black holes), and of slow 
internal motions. The main problem with this approximation is its domain of validity, which is 
limited to the near zone of the source - the region surrounding the source that is of small extent 
with respect to the wavelength of waves. A serious consequence is the a priori inability of the 
post-Newtonian expansion to incorporate the boundary conditions at infinity, which determine the 
radiation reaction force in the source's local equations of motion. The post-Minkowskian expansion, 
by contrast, is uniformly valid, as soon as the source is weakly self-gravitating, over all space-time. 
In a sense, the post-Minkowskian method is more fundamental than the post-Newtonian one; it 
can be regarded as an "upstream" approximation with respect to the post-Newtonian expansion, 

1 In this article Greek indices take the values 0, 1, 2, 3 and Latin 1, 2, 3. Our signature is +2. G and c are Newton's 
constant and the speed of light. 
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because each coefficient of the post-Minkowskian series can in turn be re-expanded in a post- 
Newtonian fashion. Therefore, a way to take into account the boundary conditions at infinity in 
the post-Newtonian series is first to perform the post-Minkowskian expansion. Notice that the 
post-Minkowskian method is also upstream (in the previous sense) with respect to the multipole 
expansion, when considered outside the source, and with respect to the far-zone expansion, when 
considered far from the source. 

The most "downstream" approximation that we shall use in this article is the post-Newtonian 
one; therefore this is the approximation that dictates the allowed physical properties of our matter 
source. We assume mainly that the source is at once slowly moving and weakly stressed, and 
we abbreviate this by saying that the source is post-Newtonian. For post-Newtonian sources, the 
parameter defined from the components of the matter stress-energy tensor T a @ and the source's 
Newtonian potential U by 



T 
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(1) 



is much less than one. This parameter represents essentially a slow motion estimate e ~ v/c, 
where v denotes a typical internal velocity. By a slight abuse of notation, following Chandrasekhar 
et al. 1681 167 ) . we shall henceforth write e = 1/c, even though e is dimensionless whereas c 
has the dimension of a velocity. The small post-Newtonian remainders will be denoted 0(l/c n ). 
Thus, 1/c <C 1 in the case of post-Newtonian sources. We have lU/c 2 ] 1 ^ 2 <C 1/c for sources 
with negligible self-gravity, and whose dynamics are therefore driven by non-gravitational forces. 
However, we shall generally assume that the source is self-gravitating; in that case we see that it 
is necessarily weakly (but not negligibly) self-gravitating, i.e. \U jc 2 \^l 2 — 0(1/ c). Note that the 
adjective "slow-motion" is a bit clumsy because we shall in fact consider very relativistic sources 
such as inspiralling compact binaries, for which 1/c can be as large as 30% in the last rotations, 
and whose description necessitates the control of high post-Newtonian approximations. 

The lowest-order wave generation formalism, in the Newtonian limit 1/c — > 0, is the famous 
quadrupole formalism of Einstein |1(J5) and Landau and Lifchitz |15Uj . This formalism can also 
be referred to as Newtonian because the evolution of the quadrupole moment of the source is 
computed using Newton's laws of gravity. It expresses the gravitational field hj^ in a transverse 
and traceless (TT) coordinate system, covering the far zone of the source 2 , as 

W - H^(N) {^(T - RIO + C (i)} + C (i) , (2, 

where R = |X| is the distance to the source, N = X/i? is the unit direction from the source to the 
observer, and V l j a b = 'Pia'Pjb - ^ $ij*Pij r Pab is the TT projection operator, with *Pij — 6%j ~- 
being the projector onto the plane orthogonal to N. The source's quadrupole moment takes the 
familiar Newtonian form 



Qij(t) = / d 3 xp(x,i) ( XiXj - -StjX 2 V 

J source \ / 



(3) 



where p is the Newtonian mass density. The total gravitational power emitted by the source in all 
directions is given by the Einstein quadrupole formula 



G f d 3 Q ab d s Q ab 
5c 5 1 dT 3 dT 3 \c 2 )j 



(4) 



2 The TT coordinate system can be extended to the near zone of the source as well; see for instance Ref. 11481 . 
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Our notation L stands for the total gravitational "luminosity" of the source. The cardinal virtues 
of the Einstein-Landau-Lifchitz quadrupole formalism are its generality - the only restrictions are 
that the source be Newtonian and bounded - its simplicity, as it necessitates only the computation 
of the time derivatives of the Newtonian quadrupole moment (using the Newtonian laws of motion), 
and, most importantly, its agreement with the observation of the dynamics of the Hulse- Taylor 
binary pulsar PSR 1913+16 |2()51 12061 1204) . Indeed the prediction of the quadrupole formalism for 
the waves emitted by the binary pulsar system comes from applying Equation (@J to a system of 
two point masses moving on an eccentric orbit (the classic reference is Peters and Mathews |175) : 
see also Refs. |108l 1213] ~). Then, relying on the energy equation 

where E is the Newtonian binary's center-of-mass energy, we deduce from Kepler's third law the 
expression of the "observable" , that is, the change in the orbital period P of the pulsar, or P, as a 
function of P itself. From the binary pulsar test, we can say that the post-Newtonian corrections to 
the quadrupole formalism, which we shall compute in this article, have already received, in the case 
of compact binaries, strong observational support (in addition to having, as we shall demonstrate, 
a sound theoretical basis). 

The multipole expansion is one of the most useful tools of physics, but its use in general relativity 
is difficult because of the non-linearity of the theory and the tensorial character of the gravitational 
interaction. In the stationary case, the multipole moments are determined by the expansion of the 
metric at spatial infinity [1171 11261 1198| . while, in the case of non-stationary fields, the moments, 
starting with the quadrupole, are defined at future null infinity. The multipole moments have 
been extensively studied in the linearized theory, which ignores the gravitational forces inside the 
source. Early studies have extended the formula (@J) to include the currcnt-quadrupolc and mass- 
octupole moments |168II167] . and obtained the corresponding formulas for linear momentum |168l 
11671 11UI 1183) and angular momentum |174l 175) . The general structure of the infinite multipole 
series in the linearized theory was investigated by several works [1881 11891 11781 12U7] , from which it 
emerged that the expansion is characterized by two and only two sets of moments: mass-type and 
current-type moments. Below we shall use a particular multipole decomposition of the linearized 
(vacuum) metric, parametrized by symmetric and trace-free (STF) mass and current moments, as 
given by Thornc 207 . The explicit expressions of the multipole moments (for instance in STF 
guise) as integrals over the source, valid in the linearized theory but irrespective of a slow motion 
hypothesis, are completely known |156l IB*4l 159] , 

In the full non-linear theory, the (radiative) multipole moments can be read off the coefficient 
of 1/R in the expansion of the metric when R — ► +oo, with a null coordinate T— R/c — const. The 
solutions of the field equations in the form of a far-field expansion (power series in 1/R) have been 
constructed, and their properties elucidated, by Bondi et al. [S3] and Sachs [19(Jj . The precise way 
under which such radiative space-times fall off asymptotically has been formulated geometrically 
by Penrose |172l 1173] in the concept of an asymptotically simple space-time (see also Ref. |118] 1. 
The resulting Bondi-Sachs-Penrose approach is very powerful, but it can answer a priori only a 
part of the problem, because it gives information on the field only in the limit where R — ► +oo, 
which cannot be connected in a direct way to the actual behaviour of the source. In particular the 
multipole moments that one considers in this approach are those measured at infinity - we call 
them the radiative multipole moments. These moments arc distinct, because of non-linearities, 
from some more natural source multipole moments, which are defined operationally by means of 
explicit integrals extending over the matter and gravitational fields. 

An alternative way of defining the multipole expansion within the complete non-linear theory is 
that of Blanchet and Damour [261112) . following pioneering work by Bonnor and collaborators 
021 OHO 1127] and Thorne )207j . In this approach the basic multipole moments are the source 
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moments, rather than the radiative ones. In a first stage, the moments are left unspecified, as being 
some arbitrary functions of time, supposed to describe an actual physical source. They are iterated 
by means of a post-Minkowskian expansion of the vacuum field equations (valid in the source's 
exterior). Technically, the post-Minkowskian approximation scheme is greatly simplified by the 
assumption of a multipolar expansion, as one can consider separately the iteration of the different 
multipole pieces composing the exterior field (whereas, the direct attack of the post-Minkowskian 
expansion, valid at once inside and outside the source, faces some calculational difficulties |212ll7fi) \ 
In this "multipolar-post-Minkowskian" formalism, which is physically valid over the entire weak- 
field region outside the source, and in particular in the wave zone (up to future null infinity), the 
radiative multipole moments are obtained in the form of some non-linear functionals of the more 
basic source moments. A priori, the method is not limited to post-Newtonian sources, however we 
shall see that, in the current situation, the closed-form expressions of the source multipole moments 
can be established only in the case where the source is post- Newtonian |15l |2"U) . The reason is that 
in this case the domain of validity of the post-Newtonian iteration (viz. the near zone) overlaps the 
exterior weak-field region, so that there exists an intermediate zone in which the post-Newtonian 
and multipolar expansions can be matched together. This is a standard application of the method 
of matched asymptotic expansions in general relativity (f>HI lf>2j . 

To be more precise, we shall show how a systematic multipolar and post-Minkowskian itera- 
tion scheme for the vacuum Einstein field equations yields the most general physically admissible 
solution of these equations 26 . The solution is specified once we give two and only two sets of 
time-varying (source) multipole moments. Some general theorems about the near-zone and far- 
zone expansions of that general solution will be stated. Notably, we find ^2 that the asymptotic 
behaviour of the solution at future null infinity is in agreement with the findings of the Bondi- 
Sachs-Penrose (531 11901 11721 11731 1118) approach to gravitational radiation. However, checking 
that the asymptotic structure of the radiative field is correct is not sufficient by itself, because 
the ultimate aim is to relate the far field to the properties of the source, and we are now obliged 
to ask: What are the multipole moments corresponding to a given stress-energy tensor T a/3 de- 
scribing the source? Only in the case of post-Newtonian sources has it been possible to answer 
this question. The general expression of the moments was obtained at the level of the second 
post-Newtonian (2PN) order in Ref. ^SJj an d was subsequently proved to be in fact valid up to 
any post-Newtonian order in Ref. |20|. The source moments are given by some integrals extending 
over the post-Newtonian expansion of the total (pseudo) stress-energy tensor r Q/3 , which is made 
of a matter part described by T a @ and a crucial non-linear gravitational source term A a @ . These 
moments carry in front a particular operation of taking the finite part (TV as we call it below) , 
which makes them mathematically well-defined despite the fact that the gravitational part A"^ 
has a spatially infinite support, which would have made the bound of the integral at spatial in- 
finity singular (of course the finite part is not added a posteriori to restore the well-definiteness 
of the integral, but is proved to be actually present in this formalism). The expressions of the 
moments had been obtained earlier at the 1PN level, albeit in different forms, in Ref. [2E1 for the 
mass-type moments (strangely enough, the mass moments admit a compact-support expression at 
1PN order), and in Ref. [!|U] f° r the current-type ones. 

The wave-generation formalism resulting from matching the exterior multipolar and post- 
Minkowskian field |26l 112] to the post-Newtonian source ^3 120] is able to take into account, 
in principle, any post-Newtonian correction to both the source and radiative multipole moments 
(for any multipolarity of the moments). The relationships between the radiative and source mo- 
ments include many non-linear multipole interactions, because the source moments mix with each 
other as they "propagate" from the source to the detector. Such multipole interactions include the 
famous effects of wave tails, corresponding to the coupling between the non-static moments with 
the total mass M of the source. The non-linear multipole interactions have been computed within 
the present wave-generation formalism up to the 3PN order in Refs. |29U21llTU] . Furthermore, the 
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back-reaction of the gravitational-wave emission onto the source, up to the 1.5PN order relative 
to the leading order of radiation reaction, has also been studied within this formalism |27U14llT%] . 
Now, recall that the leading radiation reaction force, which is quadrupolar, occurs already at the 
2.5PN order in the source's equations of motion. Therefore the 1.5PN "relative" order in the 
radiation reaction corresponds in fact to the 4PN order in the equations of motion, beyond the 
Newtonian acceleration. It has been shown that the gravitational wave tails enter the radiation re- 
action at precisely the 1.5PN relative order, which means 4PN "absolute" order [22]. A systematic 
post-Newtonian iteration scheme for the near-zone field, formally taking into account all radiation 
reaction effects, has been recently proposed, consistent with the present formalism |182l I41| . 

A different wave-generation formalism has been devised by Will and Wiseman |217) (see also 
Refs. [2HD EDI EH]), after earlier attempts by Epstein and Wagoner [HI] and Thorne [27T7] . 
This formalism has exactly the same scope as ours, i.e. it applies to any isolated post-Newtonian 
sources, but it differs in the definition of the source multipole moments and in many technical 
details when properly implemented |217) . In both formalisms, the moments are generated by the 
post-Newtonian expansion of the pseudo-tensor T a/3 , but in the Will- Wiseman formalism they 
are defined by some compact- support integrals terminating at some finite radius 1Z enclosing the 
source, e.g., the radius of the near zone). By contrast, in our case |15l 120) . the moments are 
given by some integrals covering the whole space and regularized by means of the finite part TV . 
We shall prove the complete equivalence, at the most general level, between the two formalisms. 
What is interesting about both formalisms is that the source multipole moments, which involve a 
whole series of relativistic corrections, are coupled together, in the true non-linear solution, in a 
very complicated way. These multipole couplings give rise to the many tail and related non-linear 
effects, which form an integral part of the radiative moments at infinity and thereby of the observed 
signal. 

Part [I] of this article is devoted to a presentation of the post- Newtonian wave generation for- 
malism. We try to state the main results in a form that is simple enough to be understood without 
the full details, but at the same time we outline some of the proofs when they present some interest 
on their own. To emphasize the importance of some key results, we present them in the form of 
mathematical theorems. 

1.2 Problem posed by compact binary systems 

Inspiralling compact binaries, containing neutron stars and/or black holes, are promising sources 
of gravitational waves detectable by the detectors LIGO, VIRGO, GEO and TAMA. The two 
compact objects steadily lose their orbital binding energy by emission of gravitational radiation; 
as a result, the orbital separation between them decreases, and the orbital frequency increases. 
Thus, the frequency of the gravitational-wave signal, which equals twice the orbital frequency for 
the dominant harmonics, "chirps" in time (i.e. the signal becomes higher and higher pitched) until 
the two objects collide and merge. 

The orbit of most inspiralling compact binaries can be considered to be circular, apart from the 
gradual inspiral, because the gravitational radiation reaction forces tend to circularize the motion 
rapidly. For instance, the eccentricity of the orbit of the Hulse-Taylor binary pulsar is presently 
eo = 0.617. At the time when the gravitational waves emitted by the binary system will become 
visible by the detectors, i.e. when the signal frequency reaches about 10 Hz (in a few hundred 
million years from now), the eccentricity will be e = 5.3 x 10 -6 - a value calculated from the 
Peters |174) law, which is itself based on the quadrupole formula J2Jl- 

The main point about modelling the inspiralling compact binary is that a model made of two 
structureless point particles, characterized solely by two mass parameters mi and mi (and possibly 
two spins), is sufficient. Indeed, most of the non-gravitational effects usually plaguing the dynamics 
of binary star systems, such as the effects of a magnetic field, of an interstellar medium, and so on, 
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arc dominated by gravitational effects. However, the real justification for a model of point particles 
is that the effects due to the finite size of the compact bodies are small. Consider for instance the 
influence of the Newtonian quadrupole moments Qi and Q2 induced by tidal interaction between 
two neutron stars. Let a\ and a 2 be the radius of the stars, and L the distance between the two 
centers of mass. We have, for tidal moments, 

a 5 a 5 
Qi=kim2jj, Q2 = k 2 m 1 j^, (6) 

where k\ and k 2 are the star's dimensionless (second) Love numbers |159| . which depend on their 
internal structure, and are, typically, of the order unity. On the other hand, for compact objects, 
we can introduce their "compactness" , defined by the dimensionless ratios 

Ki = -, K 2 = 5-, (7) 

a\C A a 2 & 

which equal ~ 0.2 for neutron stars (depending on their equation of state). The quadrupoles Qi and 
Q2 will affect both sides of Equation (JjJJ, i.e. the Newtonian binding energy E of the two bodies, and 
the emitted total gravitational flux £ as computed using the Newtonian quadrupole formula 
It is known that for inspiralling compact binaries the neutron stars are not co-rotating because 
the tidal synchronization time is much larger than the time left till the coalescence. As shown by 
Kochanek |144j the best models for the fluid motion inside the two neutron stars are the so-called 
Roche-Riemann ellipsoids, which have tidally locked figures (the quadrupole moments face each 
other at any instant during the inspiral) , but for which the fluid motion has zero circulation in the 
inertial frame. In the Newtonian approximation we find that within such a model (in the case of 
two identical neutron stars) the orbital phase, deduced from Equation (J5J, reads 



* finite size / 

b - 0o = 



8x 5 / 2 



1 + const k 



where x = (Gmui/c 3 ) 2 ^ 3 is a standard dimensionless post-Newtonian parameter ~ 1/c 2 (u> is the 
orbital frequency), and where k is the Love number and K is the compactness of the neutron 
star. The first term in the right-hand side of Equation (JSJ corresponds to the gravitational-wave 
damping of two point masses; the second term is the finite-size effect, which appears as a relative 
correction, proportional to (x/K) 5 , to the latter radiation damping effect. Because the finite-size 
effect is purely Newtonian, its relative correction ~ (x/K) 5 should not depend on c; and indeed 
the factors 1/c 2 cancel out in the ratio x/K. However, the compactness K of compact objects is by 
Equation (JJJ of the order unity (or, say, one half), therefore the 1/c 2 it contains should not be taken 
into account numerically in this case, and so the real order of magnitude of the relative contribution 
of the finite-size effect in Equation JSJl is given by x 5 alone. This means that for compact objects the 
finite-size effect should be comparable, numerically, to a post-Newtonian correction of magnitude 
x 5 ~ 1/c 10 namely 5PN order 3 . This is a much higher post-Newtonian order than the one at which 
we shall investigate the gravitational effects on the phasing formula. Using k' = const fc ~ 1 and 
K ~ 0.2 for neutron stars (and the bandwidth of a VIRGO detector between 10 Hz and 1000 Hz), 
we find that the cumulative phase error due to the finite-size effect amounts to less that one orbital 
rotation over a total of ~ 16, 000 produced by the gravitational-wave damping of point masses. 
The conclusion is that the finite-size effect can in general be neglected in comparison with purely 
gravitational-wave damping effects. But note that for non-compact or moderately compact objects 
(such as white dwarfs for instance) the Newtonian tidal interaction dominates over the radiation 
damping. 

3 See Ref. 1811 for the proof of such an "effacement" principle in the context of relativistic equations of motion. 
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The inspiralling compact binaries are ideally suited for application of a high-order post-Newtonian 
wave generation formalism. The main reason is that these systems are very relativistic, with 
orbital velocities as high as 0.5c in the last rotations (as compared to ~ 10~ 3 c for the binary 
pulsar) , and it is not surprising that the quadrupole-moment formalism (0 |3 21 constitutes a 
poor description of the emitted gravitational waves, since many post-Newtonian corrections play 
a substantial role. This expectation has been confirmed in recent years by several measurement- 
analyses |Z7|[ZHIIiniCnil2!S3IIHnilIHIlinSIE21, which have demonstrated that the post-Newtonian 
precision needed to implement successively the optimal filtering technique in the LIGO/ VIRGO 
detectors corresponds grossly, in the case of neutron-star binaries, to the 3PN approximation, or 
1/c 6 beyond the quadrupole moment approximation. Such a high precision is necessary because of 
the large number of orbital rotations that will be monitored in the detector's frequency bandwidth 
(~ 16,000 in the case of neutron stars), giving the possibility of measuring very accurately the 
orbital phase of the binary. Thus, the 3PN order is required mostly to compute the time evolution 
of the orbital phase, which depends, via the energy equation Q, on the center-of-mass binding 
energy E and the total gravitational-wave energy flux C. 

In summary, the theoretical problem posed by inspiralling compact binaries is two-fold: On the 
one hand E, and on the other hand £, are to be deduced from general relativity with the 3PN 
precision or better. To obtain E we must control the 3PN equations of motion of the binary in the 
case of general, not necessarily circular, orbits. As for L it necessitates the application of a 3PN 
wave generation formalism (actually, things are more complicated because the equations of motion 
are also needed during the computation of the flux) . It is quite interesting that such a high order 
approximation as the 3PN one should be needed in preparation for LIGO and VIRGO data analysis. 
As we shall see, the signal from compact binaries contains at the 3PN order the signature of several 
non-linear effects which are specific to general relativity. Therefore, we have here the possibility of 
probing, experimentally, some aspects of the non-linear structure of Einstein's theory |47l 148) . 

1.3 Post-Newtonian equations of motion and radiation 

By equations of motion we mean the explicit expression of the accelerations of the bodies in terms 
of the positions and velocities. In Newtonian gravity, writing the equations of motion for a system 
of N particles is trivial; in general relativity, even writing the equations in the case N — 2 is 
difficult. The first relativistic term, at the 1PN order, was derived by Lorentz and Droste |153) . 
Subsequently, Einstein, Infeld and Hoffmann |106j obtained the 1PN corrections by means of their 
famous "surface-integral" method, in which the equations of motion are deduced from the vacuum 
field equations, and which are therefore applicable to any compact objects (be they neutron stars, 
black holes, or, perhaps, naked singularities). The lPN-accurate equations were also obtained, for 
the motion of the centers of mass of extended bodies, by Petrova |176) and Fock |112| (see also 
Ref. [TEEp. 

The 2PN approximation was tackled by Ohta et al. |162l 11641 1163) . who considered the post- 
Newtonian iteration of the Hamiltonian of N point-particles. We refer here to the Hamiltonian 
as the Fokker-type Hamiltonian, which is obtained from the matter-plus-field Arnowitt-Deser- 
Misner (ADM) Hamiltonian by eliminating the field degrees of freedom. The result for the 2PN 
and even 2.5PN equations of binary motion in harmonic coordinates was obtained by Damour 
and Deruelle 86 85 104 80 81 , building on a non-linear iteration of the metric of two particles 
initiated in Ref. jTJ. The corresponding result for the ADM-Hamiltonian of two particles at the 
2PN order was given in Ref. [Ml (see also Refs. [H21 ESI)- Kopeikin [HEJ derived the 2.5PN 
equations of motion for two extended compact objects. The 2. 5PN- accurate harmonic-coordinate 
equations as well as the complete gravitational field (namely the metric gap) generated by two 
point masses were computed in Ref. 02] > following a method based on previous work on wave 
generation [T5| . 
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Up to the 2PN level the equations of motion are conservative. Only at the 2.5PN order ap- 
pears the first non-conservative effect, associated with the gravitational radiation reaction. The 
(harmonic-coordinate) equations of motion up to that level, as derived by Damour and Deru- 
elle jSHJ 11041 1801 181) . have been used for the study of the radiation damping of the binary 
pulsar - its orbital P |81l 1821 H()2) . It is important to realize that the 2.5PN equations of motion 
have been proved to hold in the case of binary systems of strongly self-gravitating bodies jHJ- This 
is via an "effacing" principle (in the terminology of Damour |81| ) for the internal structure of the 
bodies. As a result, the equations depend only on the "Schwarzschild" masses, mi and m^, of the 
compact objects. Notably their compactness parameters K\ and K2, defined by Equation J7J), do 
not enter the equations of motion, as has been explicitly verified up to the 2.5PN order by Kopeikin 
et al. |14fil 1124] . who made a "physical" computation, a la Fock, taking into account the internal 
structure of two self-gravitating extended bodies. The 2.5PN equations of motion have also been 
established by Itoh, Futamase and Asada [1311 1132] , who use a variant of the surface- integral ap- 
proach of Einstein, Infeld and Hoffmann |106| . that is valid for compact bodies, independently of 
the strength of the internal gravity. 

The present state of the art is the 3PN approximation 4 . To this order the equations have 
been worked out independently by two groups, by means of different methods, and with equivalent 
results. On the one hand, Jaranowski and Schafer [1361 11371 1138) . and Damour, Jaranowski, and 
Schafer [HSI EZI following the line of research of Refs. [THl E2I E31 EE] , employ the ADM- 
Hamiltonian formalism of general relativity; on the other hand, Blanchet and Faye |37l 1381 151)1 139| . 
and de Andrade, Blanchet, and Faye |103) . founding their approach on the post-Newtonian itera- 
tion initiated in Ref. [42] . compute directly the equations of motion (instead of a Hamiltonian) in 
harmonic coordinates. The end results have been shown [9711103) to be physically equivalent in the 
sense that there exists a unique "contact" transformation of the dynamical variables that changes 
the harmonic-coordinates Lagrangian obtained in Ref. [103] into a new Lagrangian, whose Legen- 
dre transform coincides exactly with the Hamiltonian given in Ref. |95j . The 3PN equations of 
motion, however, depend on one unspecified numerical coefficient, w s t a tic in the ADM-Hamiltonian 
formalism and A in the harmonic-coordinates approach, which is due to some incompleteness of 
the Hadamard self-field regularization method. This coefficient has been fixed by means of a 
dimensional regularization, both within the ADM-Hamiltonian formalism |96j . and the harmonic- 
coordinates equations of motion |30) . The works |96l 130) have demonstrated the power of dimen- 
sional regularization and its perfect adequateness for the problem of the interaction between point 
masses in general relativity. Furthermore, an important work by Itoh and Futamase |130II129] (us- 
ing the same surface-integral method as in Refs. [1311 1132] ) succeeded in obtaining the complete 
3PN equations of motion in harmonic coordinates directly, i.e. without ambiguity and containing 
the correct value for the parameter A. 

So far the status of the post-Newtonian equations of motion is quite satisfying. There is mutual 
agreement between all the results obtained by means of different approaches and techniques, when- 
ever it is possible to compare them: point particles described by Dirac delta- functions, extended 
post-Newtonian fluids, surface-integrals methods, mixed post-Minkowskian and post-Newtonian 
expansions, direct post-Newtonian iteration and matching, harmonic coordinates versus ADM- 
type coordinates, and different processes or variants of the regularization of the self field of point 
particles. In Part [H] of this article, we shall present the complete results for the 3PN equations of 
motion, and for the associated Lagrangian and Hamiltonian formulations (from which we deduce 
the center-of-mass energy E). 

The second sub-problem, that of the computation of the energy flux £, has been carried out by 

4 Let us mention that the 3.5PN terms in the equations of motion are also known, both for point-particle bina- 
ries IT^HTMIIT^5HT7THTI5HTrTTl and extended fluid bodies Q^IEll; they correspond to 1PN "relative" corrections 
in the radiation reaction force. Known also is the contribution of wave tails in the equations of motion, which arises 
at the 4PN order and represents a 1.5PN modification of the gravitational radiation damping |27| . 
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application of the wave-generation formalism described previously. Following earliest computations 
at the 1PN level |214l I49j . at a time when the post-Newtonian corrections in C had a purely 
academic interest, the energy flux of inspiralling compact binaries was completed to the 2PN order 
by Blanchet, Damour and Iyer |33l 1119) . and, independently, by Will and Wiseman |217) . using 
their own formalism (see Refs. EEI for joint reports of these calculations). The preceding 
approximation, 1.5PN, which represents in fact the dominant contribution of tails in the wave 
zone, had been obtained in Refs. [2181 ED] by application of the formula for tail integrals given in 
Ref. [2H|- Higher-order tail effects at the 2.5PN and 3.5PN orders, as well as a crucial contribution 
of tails generated by the tails themselves (the so-called "tails of tails") at the 3PN order, were 
obtained by Blanchet [TH1 CSJ. However, unlike the 1.5PN, 2.5PN, and 3.5PN orders that are 
entirely composed of tail terms, the 3PN approximation also involves, besides the tails of tails, 
many non-tail contributions coming from the rclativistic corrections in the (source) multipole 
moments of the binary. These have been "almost" completed in Refs. [451 14UI |4"4"] . in the sense 
that the result still involves one unknown numerical coefficient, due to the use of the Hadamard 
regularization, which is a combination of the parameter A in the equations of motion, and a new 
parameter 9 coming from the computation of the 3PN quadrupole moment. The latter parameter 
is itself a linear combination of three unknown parameters, 9 = £ + 2k + £. We shall review the 
computation of the three parameters £, k, and £ by means of dimensional regularization [311 132) . 
In Part[H]of this article, we shall present the most up-to-date results for the 3.5PN energy flux and 
orbital phase, deduced from the energy balance equation J5J, supposed to be valid at this order. 

The post-Newtonian flux C, which comes from a "standard" post-Newtonian calculation, is in 
complete agreement (up to the 3.5PN order) with the result given by the very different technique 
of linear black-hole perturbations, valid in the "test-mass" limit where the mass of one of the 
bodies tends to zero (limit v — > 0, where v = fi/m). Linear black-hole perturbations, triggered by 
the geodesic motion of a small mass around the black hole, have been applied to this problem by 
Poisson 1 179] at the 1.5PN order (following the pioneering work of Galt'sov et al. |llfij ). and by 
Tagoshi and Nakamura |2()0) . using a numerical code, up to the 4PN order. This technique has 
culminated with the beautiful analytical methods of Sasaki, Tagoshi and Tanaka [1911 12021 12()3| 
(see also Ref. [157] 1. who solved the problem up to the extremely high 5.5PN order. 
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Part I 

Post-Newtonian Sources 



2 Einstein's Field Equations 

The field equations of general relativity form a system of ten second-order partial differential 
equations obeyed by the space-time metric g a j3, 

G* > \g,dg,&g] = ^T#\g], (9) 

where the Einstein curvature tensor G Q/3 = i? Q/3 — ^Rg a @ is generated, through the gravitational 
coupling k — 8irG/c 4 , by the matter stress-energy tensor T™^. Among these ten equations, four 
govern, via the contracted Bianchi identity, the evolution of the matter system, 

V M G Q " = => V A1 T^ = 0. (10) 

The space-time geometry is constrained by the six remaining equations, which place six independent 
constraints on the ten components of the metric g a p, leaving four of them to be fixed by a choice 
of a coordinate system. 

In most of this paper we adopt the conditions of harmonic, or de Donder, coordinates. We 
define, as a basic variable, the gravitational-field amplitude 

h aP = V=gg a0 - v aP , (ii) 

where g a @ denotes the contravariant metric (satisfying g a>l gp.fi = 5g), where g is the determinant of 
the covariant metric, g — det(<? Q ^), and where rj a ^ represents an auxiliary Minkowskian metric. The 
harmonic-coordinate condition, which accounts exactly for the four equations Ijl0|l corresponding 
to the conservation of the matter tensor, reads 

d^h al1 = 0. (12) 

Equations I12|) introduce into the definition of our coordinate system a preferred Minkowskian 
structure, with Minkowski metric rj a p. Of course, this is not contrary to the spirit of general rel- 
ativity, where there is only one physical metric g a p without any flat prior geometry, because the 
coordinates are not governed by geometry (so to speak) , but rather are chosen by researchers when 
studying physical phenomena and doing experiments. Actually, the coordinate condition (|12|) is 
especially useful when we view the gravitational waves as perturbations of space-time propagating 
on the fixed Minkowskian manifold with the background metric rj a p. This view is perfectly legiti- 
mate and represents a fruitful and rigorous way to think of the problem when using approximation 
methods. Indeed, the metric r] a p, originally introduced in the coordinate condition l|12|l . does exist 
at any finite order of approximation (neglecting higher-order terms), and plays in a sense the role 
of some "prior" flat geometry. 

The Einstein field equations in harmonic coordinates can be written in the form of inhomoge- 
ncous flat d'Alembertian equations, 

nh af > = (13) 

c 4 

where □ = □,,= if^d^dv. The source term r a/3 can rightly be interpreted as the stress-energy 
pseudo-tensor (actually, T a/3 is a Lorentz tensor) of the matter fields, described by T Q/3 , and the 
gravitational field, given by the gravitational source term A a @ , i.e. 

r aP = \g\T# + ^A^. (14) 
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The exact expression of A" 13 , including all non-linearities, reads 5 

-g^g.rdxh^d.h^ - g^g VT d x h aT d^ x + g^g XT d x h^d T h^ 

+ l -{2g^g^ - g^gn&gxrg^ ~ g rt gx«)d^d„h^. (15) 

As is clear from this expression, A a @ is made of terms at least quadratic in the gravitational-field 
strength h and its first and second space-time derivatives. In the following, for the highest post- 
Newtonian order that we consider (3PN), we need the quadratic, cubic and quartic pieces of A a °. 
With obvious notation, we can write them as 

A af3 = N a(3 [h, h] + M a/3 [h, h, h] + L al3 [h, h, h, h] + 0(h 5 ). (16) 

These various terms can be straightforwardly computed from Equation 1)15(1: see Equations (3.8) 
in Ref. |2H| for explicit expressions. 

As said above, the condition l|12|l is equivalent to the matter equations of motion, in the sense 
of the conservation of the total pseudo-tensor r"' 5 , 

=0 V^r" = 0. (17) 

In this article, we look for the solutions of the field equations H13lll4llT51ll7fl under the following 
four hypotheses: 

1. The matter stress-energy tensor T a ° is of spatially compact support, i.e. can be enclosed 
into some time-like world tube, say r < a, where r = |x| is the harmonic-coordinate radial 
distance. Outside the domain of the source, when r > a, the gravitational source term, 
according to Equation ifTTJl. is divergence- free, 

d^A afl — (whenr>a). (18) 



2. The matter distribution inside the source is smooth 6 : T Q/3 £ C°°(IR 3 ). We have in mind a 
smooth hydrodynamical "fluid" system, without any singularities nor shocks (a priori), that 
is described by some Eulerian equations including high relativistic corrections. In particular, 
we exclude from the start any black holes (however we shall return to this question when we 
find a model for describing compact objects). 

3. The source is post- Newtonian in the sense of the existence of the small parameter defined 
by Equation For such a source we assume the legitimacy of the method of matched 
asymptotic expansions for identifying the inner post-Newtonian field and the outer multipolar 
decomposition in the source's exterior near zone. 

4. The gravitational field has been independent of time (stationary) in some remote past, i.e. 
before some finite instant — T in the past, in the sense that 

[/i a/3 (x,<)] = wheni<-T. (19) 

5 See also Equation 11401 for the expression in d + 1 space-time dimensions. 

6 N, Z, R, and C are the usual sets of non-negative integers, integers, real numbers, and complex numbers; C P (Q) 
is the set of p-times continuously differentiable functions on the open domain Q (p < +00). 
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The latter condition is a means to impose, by brute force, the famous no-incoming radiation 
condition, ensuring that the matter source is isolated from the rest of the Universe and does not 
receive any radiation from infinity. Ideally, the no-incoming radiation condition should be imposed 
at past null infinity. We shall later argue (see Section |HJ that our condition of stationarity in the 
past, Equation Ijl9|l . although much weaker than the real no- incoming radiation condition, does 
not entail any physical restriction on the general validity of the formulas we derive. 

Subject to the condition p9|l. the Einstein differential field equations (|13JI can be written 
equivalently into the form of the integro-differential equations 

= ^D-lr^, (20) 
containing the usual retarded inverse d'Alembertian operator, given by 

(□-i/)(x,i)=-i- JJJ ^^/( x ',*-|x-x'|/c), (21) 

extending over the whole three-dimensional space R 3 . 
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3 Linearized Vacuum Equations 



In what follows we solve the field equations ((121 1 1 3|> . in the vacuum region outside the compact- 
support source, in the form of a formal non-linearity or post-Minkowskian expansion, considering 
the field variable h a/3 as a non-linear metric perturbation of Minkowski space-time. At the lin- 
earized level (or first-post-Minkowskian approximation), we write: 

Kg = Ghf + G(G 2 ), (22) 

where the subscript "ext" reminds us that the solution is valid only in the exterior of the source, and 
where we have introduced Newton's constant G as a book-keeping parameter, enabling one to label 
very conveniently the successive post-Minkowskian approximations. Since h al3 is a dimensionless 
variable, with our convention the linear coefficient in Equation 1(22(1 has the dimension of the 
inverse of G - a mass squared in a system of units where K = c = 1. In vacuum, the harmonic- 
coordinate metric coefficient satisfies 

Uhf = 0, (23) 
W = 0. (24) 

We want to solve those equations by means of an infinite multipolar series valid outside a time- 
like world tube containing the source. Indeed the multipole expansion is the correct method for 
describing the physics of the source as seen from its exterior (r > a). On the other hand, the post- 
Minkowskian series is physically valid in the weak-field region, which surely includes the exterior of 
any source, starting at a sufficiently large distance. For post-Newtonian sources the exterior weak- 
field region, where both multipole and post-Minkowskian expansions are valid, simply coincides 
with the exterior r > a. It is therefore quite natural, and even, one would say inescapable when 
considering general sources, to combine the post-Minkowskian approximation with the multipole 
decomposition. This is the original idea of the "double-expansion" series of Bonnor [Hi], which 
combines the G-expansion (or m-expansion in his notation) with the a-expansion (equivalent to the 
multipole expansion, since the Ith order multipole moment scales like a 1 with the source radius). 

The multipolar-post-Minkowskian method will be implemented systematically, using STF- 
harmonics to describe the multipole expansion |2()7) . and looking for a definite algorithm for the 
approximation scheme [2| . The solution of the system of equations 1(231 124(1 takes the form of a 
series of retarded multipolar waves 7 

*-_g ft («ZSL^), (25) 

where r = |x|, and where the functions K^f = K^f i are smooth functions of the retarded 
time u = t — r/c [K^(u) € C°°(R)], which become constant in the past, when t < —T. It 
is evident, since a monopolar wave satisfies D(K^(u)/r) = and the d'Alembertian commutes 

7 Our notation is the following: L = i\ii ■ ■ .i\ denotes a multi-index, made of I (spatial) indices. Similarly we 
write for instance P = ji . . . j p (in practice, we generally do not need to consider the carrier letter t or j), or 
aL — 1 = aii ■ ■ ■ H— i- Always understood in expressions such as Equation 1251 are I summations over the I indices 
ranging from 1 to 3. The derivative operator 8l is a short- hand for di 1 ...8^. The function Kl is 
symmetric and trace-free (STF) with respect to the I indices composing L. This means that for any pair of indices 
i P ,iq € L, we have JsT...i p ...i s ... = K,,.i q ,,.i and that 5i„i„K,,.i p ..,i q ... = (see Ref. 12071 and Appendices A and 
B in Ref. 26 for reviews about the STF formalism). The STF projection is denoted with a hat, so = K^, 
or sometimes with carets around the indices, Kl = K(L)- in particular, = is the STF projection of the 

product of unit vectors = rii 1 . . . ; an expansion into STF tensors % = h^ij), (f>) is equivalent to the usual 
expansion in spherical harmonics Y; m = Y; m (0,0). Similarly, we denote = Xi 1 . . . Xi l = r'n^ and = x^ L y 
Superscripts like (p) indicate p successive time-derivations. 
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with the multi-derivative <9l, that Equation 1)25|) represents the most general solution of the wave 
equation (|23|) (see Section|21in Ref. 26 for a proof based on the Euler-Poisson-Darboux equation). 
The gauge condition (|24|) . however, is not fulfilled in general, and to satisfy it we must algebraically 
decompose the set of functions K^ 1 , K l l into ten tensors which are STF with respect to all 
their indices, including the spatial indices i, ij. Imposing the condition (|24|) reduces the number 
of independent tensors to six, and we find that the solution takes an especially simple "canonical" 
form, parametrized by only two moments, plus some arbitrary linearized gauge transformation |2()7I 
Ell- 
Theorem 1 The most general solution of the linearized field equations \2tft \24[ ), outside some 
time-like world tube enclosing the source (r > a), and stationary in the past (see Equation \19t l). 
reads 

hf = kf + d a ^ + dV-^ a/3 «. (26) 

The first term depends on two STF-tensorial multipole moments, Il(u) and 3l{u), which are 
arbitrary functions of time except for the laws of conservation of the monopole: I = const, and 
dipoles: li — const, J,; = const. It is given by 




The other terms represent a linearized gauge transformation, with gauge vector ip" of the type 
and parametrized for four other multipole moments, say Wi(w), Xi(w), Yl(u) and Zl(u). 

The conservation of the lowest-order moments gives the constancy of the total mass of the source, 
M = I = const, center-of-mass position 8 , Xj = I,/I = const, total linear momentum Pj = If' = 0, 
and total angular momentum, Si = Ji = const. It is always possible to achieve X^ = by 
translating the origin of our coordinates to the center of mass. The total mass M is the ADM 
mass of the Hamiltonian formulation of general relativity. Note that the quantities M, X^, P^ and 
Si include the contributions due to the waves emitted by the source. They describe the "initial" 
state of the source, before the emission of gravitational radiation. 

The multipole functions 1l(u) and Ji(tt), which thoroughly encode the physical properties of 
the source at the linearized level (because the other moments Wz,, . . . ,Zl parametrize a gauge 
transformation), will be referred to as the mass-type and current-type source multipole moments. 
Beware, however, that at this stage the moments are not specified in terms of the stress-energy 
tensor T af3 of the source: the above theorem follows merely from the algebraic and differential 
properties of the vacuum equations outside the source. 

8 The constancy of the center of mass - rather than a linear variation with time — results from our assumption 
of stationarity before the date — T, Hence, = 0. 
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For completeness, let us give the components of the gauge- vector ipf entering Equation (|26(l : 



z>o ' v 7 

^ = -^E^7T^(^w) (28) 
z>o ' v 7 

Because the theory is covariant with respect to non-linear diffeomorphisms and not merely with 
respect to linear gauge transformations, the moments Wl, . . . , Zl do play a physical role start- 
ing at the non-linear level, in the following sense. If one takes these moments equal to zero and 
continues the calculations one ends up with a metric depending on II and Jl only, but that 
metric will not describe the same physical source as the one constructed from the six moments 
Ii, . . .,Zi. In other words, the two non- linear metrics associated with the sets of multipole mo- 
ments {II, Jl, 0, . . . , 0} and {l L , 3 L , Wl, ■ ■ ■ , Z L } are not isometric. We point out in Section ET31 
below that the full set of moments {II, Jl, Wl, ■ ■ ■ ,Zl} is in fact physically equivalent to some 
reduced set {Ml,Sl,0, ... ,0}, but with some moments Ml, Sl that differ from I L , Jl by non- 
linear corrections (see Equation HiJfif) ). All the multipole moments II, Jl, Wl, Xl, Yl, Zl will be 
computed in Section [5] 
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4 Non-linear Iteration of the Field Equations 



By Theorem ^ we know the most general solution of the linearized equations in the exterior of 
the source. We then tackle the problem of the post-Minkowskian iteration of that solution. We 
consider the full post-Minkowskian series 

+ OC 

Kg = ^G n hf, (29) 

71 = 1 

where the first term is composed of the result given by Equations 1(261 1271 l2"%f . In this article, we 
shall always understand the infinite sums such as the one in Equation l(29[l in the sense of formal 
power series, i.e. as an ordered collection of coefficients, e.g., (ft." /3 ) ngN . We do not attempt to 
control the mathematical nature of the series and refer to the mathematical-physics literature for 
discussion (in the present context, see Refs. [T^ITUmiTgllT^ITHS] ). 



4.1 The post-Minkowskian solution 

We insert the ansatz l)29|) into the vacuum Einstein field equations ((121 113(1 . i.e. with r Q/3 = 
c 4 /(WirG)A a ° , and we equate term by term the factors of the successive powers of our book- 
keeping parameter G. We get an infinite set of equations for each of the h^'s: Vn > 2, 

Dhf =Af[h!,h 2 ,..., hn-i], (30) 

W = 0. (31) 

The right-hand side of the wave equation 1(30(1 is obtained from inserting the previous iterations, up 
to the order n — 1, into the gravitational source term. In more details, the series of equations ((30|) 
reads 

Dhf = N a P[h u hi], (32) 
□/i^ = M a P[hi, hi, hi] + N a P[h\, h 2 ] + N a ^[h 2 , hi], (33) 
nh% p = L af3 [hi,hi,hi,hi] 

+ M al3 [hi,hi, h 2 ] + M a/3 [hi,h 2 ,hi] + M a/3 [h 2 ,hi,hi] 

+ N a P [h 2 ,h 2 ]+N a > 3 [h 1 , h a ] + N a ? [h a , hi] 

(34) 

The quadratic, cubic and quartic pieces of A a/3 are defined by Equation l|16(l . 

Let us now proceed by induction. Some n being given, we assume that we succeeded in 
constructing, from the linearized coefficient hi, the sequence of post-Minkowskian coefficients 
h 2 , ha, ■ ■ ■ , h n — i, and from this we want to infer the next coefficient h n . The right-hand side 
of Equation l(30j) . A"' 3 , is known by induction hypothesis. Thus the problem is that of solving a 
wave equation whose source is given. The point is that this wave equation, instead of being valid 
everywhere in R 3 , is correct only outside the matter (r > a), and it makes no sense to solve it by 
means of the usual retarded integral. Technically speaking, the right-hand side of Equation (|30(l 
is composed of the product of many multipole expansions, which are singular at the origin of the 
spatial coordinates r = 0, and which make the retarded integral divergent at that point. This does 
not mean that there are no solutions to the wave equation, but simply that the retarded integral 
does not constitute the appropriate solution in that context. 

What we need is a solution which takes the same structure as the source term A"' 3 , i.e. is 
expanded into multipole contributions, with a singularity at r — 0, and satisfies the d'Alembcrtian 
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equation as soon as r > 0. Such a particular solution can be obtained, following the suggestion 
in Ref. |26|. by means of a mathematical trick in which one first "regularizes" the source term 
K^f by multiplying it by the factor r B , where B £ C. Let us assume, for definiteness, that A"^ 
is composed of multipolar pieces with maximal multipolarity / max - This means that we start the 
iteration from the linearized metric (|26l 1271 128|l in which the multipolar sums are actually finite 9 . 
The divergences when r — > of the source term are typically power- like, say l/r k (there are also 
powers of the logarithm of r), and with the previous assumption there will exist a maximal order 
of divergency, say /c max - Thus, when the real part of B is large enough, i.e. S? (B) > k max — 3, the 
"regularized" source term r s A"^ is regular enough when r — ► so that one can perfectly apply 
the retarded integral operator. This defines the -B-dependent retarded integral 

I^(B) = n-\ [r B Af] , (35) 

where the symbol D~J stands for the retarded integral l|21|l . It is convenient to introduce inside 
the regularizing factor some arbitrary constant length scale ro in order to make it dimensionless. 
Everywhere in this article we pose 

r= — . (36) 

The fate of the constant ro in a detailed calculation will be interesting to follow, as we shall 
see, because it provides some check that the calculation is going well. Now the point for our 
purpose is that the function I a/3 (B) on the complex plane, which was originally defined only when 
•Ji (B) > /c max — 3, admits a unique analytic continuation to all values of B G C except at some 
integer values. Furthermore, the analytic continuation of I a @(B) can be expanded, when B — > 
(namely the limit of interest to us) into a Laurent expansion involving in general some multiple 
poles. The key idea, as we shall prove, is that the finite part, or the coefficient of the zeroth power 
of B in that expansion, represents the particular solution we are looking for. We write the Laurent 
expansion of I al3 (B), when B — > 0, in the form 

r*P(B) = ^2 ifB p , (37) 

P=P0 

where pgZ, and the various coefficients are functions of the field point (x, t). When po < —1 
there are poles; —po, which depends on n, refers to the maximal order of the poles. By applying 
the box operator onto both sides of Equation 1|37[) . and equating the different powers of B, we 
arrive at 

Po<P<-l => Uuf = Q, 

p >0 D^-SSdV <38) 

F p\ 

As we see, the case p — shows that the finite-part coefficient in Equation i|37f) . namely , is a 
particular solution of the requested equation: Ol^ = A" 73 . Furthermore, we can prove that this 
term, by its very construction, owns the same structure made of a multipolar expansion singular 
at r = 0. 

Let us forget about the intermediate name L^f \ and denote, from now on, the latter solution 
by u^P = , or, in more explicit terms, 

<[= TV B =* n£ [r B A^] , (39) 

9 This assumption is justified because we are ultimately interested in the radiation field at some given finite 
post-Newtonian precision like 3PN, and because only a finite number of multipole moments can contribute at any 
finite order of approximation. With a finite number of multipoles in the linearized metric 1261 1271 1281 . there is a 
maximal multipolarity Z max (n) at any post-Minkowskian order n, which grows linearly with n. 
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where the finite-part symbol TVb=§ means the previously detailed operations of considering the 
analytic continuation, taking the Laurent expansion, and picking up the finite-part coefficient when 
B — > 0. The story is not complete, however, because u" 13 does not fulfill the constraint of harmonic 
coordinates l|31[l : its divergence, say w" = c^it"^, is different from zero in general. From the fact 
that the source term is divergence-free in vacuum, d^A"^ = (see Equation ((T%)l L we find instead 



Br —K c 

r r 



(40) 



The factor B comes from the differentiation of the regularization factor r B . So, u>" is zero only 
in the special case where the Laurent expansion of the retarded integral in Equation (|40(l does not 
develop any simple pole when B — > 0. Fortunately, when it does, the structure of the pole is quite 
easy to control. We find that it necessarily consists of a solution of the source-free d'Alembertian 
equation, and, what is more (from its stationarity in the past), the solution is a retarded one. Hence, 
taking into account the index structure of u>" , there must exist four STF-tensorial functions of the 
retarded time u = t — r/c, say Nl(u), Pl(u), Ql(u) and Rl(u), such that 
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d L [r- l N L [u)\ , 



d lL [r^PLiu)] [r^QiL-M] + e iab d aL ^ [r- 1 R bL -i(u)] } 

i=i 



(41) 



From that expression we are able to find a new object, say , which takes the same structure as 
(a retarded solution of the source-free wave equation) and, furthermore, whose divergence is 
exactly the opposite of the divergence of , i.e. d^v^ — —w%. Such a v^f is not unique, but we 
shall see that it is simply necessary to make a choice for v^f (the simplest one) in order to obtain 
the general solution. The formulas that we adopt are 



d a 



1 (-Q\- 1] +3Pf >) -s iab d a [r-'R^] [r^Na-i] , 



1=2 



+ 00 s 

< = Sijr^P + J2\ 2««a £ -i [t-'Pl-i] - 69 L _ 2(l [r- 1 P, )L _ 2 ] 

1=2 ^ 



(42) 



3^1-2 - QijL-2)] - 2d aL -2 [r- l e ab{i R 



■j)bL-2\ 



Notice the presence of anti-derivatives, denoted, e.g., by N^^(u) = J*™ dvN(v); there is no 
problem with the limit v — > — oo since all the corresponding functions are zero when t < —T. The 
choice made in Equations (|42(l is dictated by the fact that the 00 component involves only some 
monopolar and dipolar terms, and that the spatial trace ii is monopolar: v™ — — 3r~ 1 P. Finally, 
if we pose 



hf=uf+vf, (43) 

we see that we solve at once the d'Alembertian equation (|30|) and the coordinate condition l|31|l . 
That is, we have succeeded in finding a solution of the field equations at the nth post-Minkowskian 
order. By induction the same method applies to any order n, and, therefore, we have constructed a 
complete post-Minkowskian series l|29(l based on the linearized approximation h" given by Equa- 
tions 1261 1271 l28fl . The previous procedure constitutes an algorithm, which could be implemented 
by an algebraic computer programme. 
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4.2 Generality of the solution 

We have a solution, but is that a general solution? The answer, yes, is provided by the following 
result p)] : 

Theorem 2 The most general solution of the harmonic- coordinates Einstein field equations in 
the vacuum region outside an isolated source, admitting some post-Minkowskian and multipolar ex- 
pansions, is given by the previous construction as h af3 — X)n=i ^"^n^P-kj ^l, ■ ■ ■ , ZjJ. It depends 
on two sets of arbitrary STF-tensorial functions of time Il(u) and 3l(u) (satisfying the conser- 
vation laws) defined by Equations \2 r /\j , and on four supplementary functions Wl(u), . . . ,Zl(u) 
parametrizing the gauge vector ]2ty). 

The proof is quite easy. With Equation l|43[) we obtained a particular solution of the system of 
equations l|3UII31|l . To it we should add the most general solution of the corresponding homogeneous 
system of equations, which is obtained by setting K^f = into Equations (|30l I31fl . But this 
homogeneous system of equations is nothing but the linearized vacuum field equations IL'-il I24|) , 
for which we know the most general solution h°[ given by Equations (|26l 1271 12"51) . Thus, we must 
add to our "particular" solution h^f a general homogeneous solution that is necessarily of the 
type . . . ,5Zl], where 5Il, . . . , 6Zl denote some "corrections" to the multipole moments 

at the nth post-Minkowskian order. It is then clear, since precisely the linearized metric is a linear 
functional of all these moments, that the previous corrections to the moments can be absorbed 
into a re-definition of the original ones II, . . . , Zl by posing 

I ncw = Ii + G n-l (5lL; (44) 

Z L ow = Z L + G n - 1 5Z L . (45) 

After re-arranging the metric in terms of these new moments, taking into account the fact that the 
precision of the metric is limited to the nth post-Minkowskian order, and dropping the superscript 
"new" , we find exactly the same solution as the one we had before (indeed, the moments are 
arbitrary functions of time) - hence the proof. 

The six sets of multipole moments Il(u), . . . , Zl(u) contain the physical information about any 
isolated source as seen in its exterior. However, as we now discuss, it is always possible to find two, 
and only two, sets of multipole moments, M^(u) and Sl(u), for parametrizing the most general 
isolated source as well. The route for constructing such a general solution is to get rid of the 
moments Wi, Xl, Yl, Zl at the linearized level by performing the linearized gauge transformation 
Sx a = ipi, where ipf is the gauge vector given by Equations Ij28(l . So, at the linearized level, 
we have only the two types of moments Mi and Si, parametrizing k"^ by the same formulas 
as in Equations JUJl. We must be careful to denote these moments with some names different 
from 1^ and Jl because they will ultimately correspond to a different physical source. Then 
we apply exactly the same post-Minkowskian algorithm, following the formulas (|39l 1401 1411 1421 
I43fl as we did above, but starting from the gauge-transformed linear metric k^ instead of h"^ . 
The result of the iteration is therefore some k a " = E,triG"^[M L ,S L ]. Obviously this post- 
Minkowskian algorithm yields some simpler calculations as we have only two multipole moments 
to iterate. The point is that one can show that the resulting metric fc^jM^, Sl] is isometric to the 
original one Jl, • ■ ■ , Zl] if and only if M^ and Sl are related to the moments 1l, Ji, • ■ • , Zl 

by some (quite involved) non- linear equations. Therefore, the most general solution of the field 
equations, modulo a coordinate transformation, can be obtained by starting from the linearized 
metric ^"^[MljSl] instead of the more complicated fcf^pL, Jl] + d a <fi + d^ip" — i^^d^ipi, and 
continuing the post-Minkowskian calculation. 
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So why not consider from the start that the best description of the isolated source is provided 
by only the two types of multipole moments, Mi and Si, instead of the six, Ii, Ji, . . . , Zi? The 
reason is that we shall determine (in Theorem below) the explicit closed-form expressions of 
the six moments Ii, Ji, . . . , Zi, but that, by contrast, it seems to be impossible to obtain some 
similar closed- form expressions for Ml and Si. The only thing we can do is to write down the 
explicit non-linear algorithm that computes M^, Si starting from Ii, Ji, . . . , Zi. In consequence, 
it is better to view the moments Ii, Jl, • ■ ■ > Zl as more "fundamental" than Mi and Si, in the 
sense that they appear to be more tightly related to the description of the source, since they 
admit closed-form expressions as some explicit integrals over the source. Hence, we choose to refer 
collectively to the six moments Ii, Ji, . . . , Zi as the multipole moments of the source. This being 
said, the moments Mi and Si are often useful in practical computations because they yield a 
simpler post-Minkowskian iteration. Then, one can generally come back to the more fundamental 
source-rooted moments by using the fact that Mi and Si differ from the corresponding Ii and 
Ji only by high-order post-Newtonian terms like 2.5PN; see Ref. |flj)J and Equation l)96[l below. 
Indeed, this is to be expected because the physical difference between both types of moments stems 
only from non-linearities. 

4.3 Near-zone and far-zone structures 

In our presentation of the post-Minkowskian algorithm ([391 1401 1411 1421 143ft we have omitted a 
crucial recursive hypothesis, which is required in order to prove that at each post-Minkowskian 
order n, the inverse d'Alembertian operator can be applied in the way we did (and notably that the 
independent retarded integral can be analytically continued down to a neighbourhood of B = 0) . 
This hypothesis is that the "near-zone" expansion, i.e. when r — > 0, of each one of the post- 
Minkowskian coefficients h%f has a certain structure. This hypothesis is established as a theorem 
once the mathematical induction succeeds. 

Theorem 3 The general structure of the expansion of the post-Minkowskian exterior metric in 
the near-zone (when r — > 0J is of the type: VA G N, 

Mx, *) = n L r m {\nrYF L>m , v , n {t) + o(r N ), (46) 

where m £ Z, with mo < m < N (and mo becoming more and more negative as n grows), peN 
with p < n— 1. The functions F^^ rn , v , n are multilinear functionals of the source multipole moments 
Ii, . . . , Zi. 

For the proof see Ref. [22] 10 . As we see, the near-zone expansion involves, besides the simple 
powers of r, some powers of the logarithm of r, with a maximal power of n — 1. As a corollary of 
that theorem, we find (by restoring all the powers of c in Equation l[46[i and using the fact that 
each r goes into the combination r/c), that the general structure of the post-Newtonian expansion 
(c — * +oo) is necessarily of the type 

p,gGN 

where p < n — 1 (and q > 2). The post-Newtonian expansion proceeds not only with the normal 
powers of 1/c but also with powers of the logarithm of c [26) . 

Paralleling the structure of the near-zone expansion, we have a similar result concerning the 
structure of the far-zone expansion at Minkowskian future null infinity, i.e. when r — » +oo with 

10 The o and O Landau symbols for remainders have their standard meaning. 
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u = t - r/c = const: VW € N, 



h n (x, i) = J] " L ^ r)P GL, t|P ,n(«) + o f -Lj , (48) 

where fc,p S N, with 1 < k < JV, and where, likewise in the near-zone expansion l|46|l . some 
powers of logarithms, such that p < n — 1, appear. The appearance of logarithms in the far-zone 
expansion of the harmonic-coordinates metric has been known since the work of Fock |113) . One 
knows also that this is a coordinate effect, because the study of the "asymptotic" structure of 
space-time at future null infinity by Bondi et al. Sachs |19U) . and Penrose [1721 1173*] . has 
revealed the existence of other coordinate systems that avoid the appearance of any logarithms: 
the so-called radiative coordinates, in which the far-zone expansion of the metric proceeds with 
simple powers of the inverse radial distance. Hence, the logarithms are simply an artifact of the 
use of harmonic coordinates |128l 1154] , The following theorem, proved in Ref. J2]> shows that 
our general construction of the metric in the exterior of the source, when developed at future null 
infinity, is consistent with the Bondi-Sachs-Penrose [SHI I19UI 11721 1173) approach to gravitational 
radiation. 

Theorem 4 The most general multipolar-post-Minkowskian solution, stationary in the past (see 
Equation f 19$), admits some radiative coordinates (T, X), for which the expansion at future null 
infinity, R — ► +oo with U = T — R/c = const, takes the form 

H n (X, T) = J2 ^K L , k , n (U) + O (Jp) . (49) 

The functions Kz,^,n computable functionals of the source multipole moments. In radiative 
coordinates the retarded time U = T — R/c is a null coordinate in the asymptotic limit. The metric 
Hext = En>i G n H"P is asymptotically simple in the sense of Penrose \1 7lfl [TT^il , perturbatively to 
any post-Minkowskian order. 

Proof: We introduce a linearized "radiative" metric by performing a gauge transformation of the 
harmonic-coordinates metric defined by Equations (|26l 1271 128|) , namely 

Hf = hf + d a ^ + - v afi d^i, (50) 

where the gauge vector £f is 

= 2M77 OQ ln (^] . (51) 
This gauge transformation is non-harmonic: 



ro 



dpH?" = = ™ V 0a . (52) 

Its effect is to "correct" for the well-known logarithmic deviation of the harmonic coordinates' 
retarded time with respect to the true space-time characteristic or light cones. After the change of 
gauge, the coordinate u = t — r/c coincides with a null coordinate at the linearized level 11 . This 
is the requirement to be satisfied by a linearized metric so that it can constitute the linearized 
approximation to a full (post-Minkowskian) radiative field |154) . One can easily show that, at the 
dominant order when r — > +oo, 

k^KH^ =o(\\ (53) 



11 In this proof the coordinates are considered as dummy variables denoted (t, r). At the end, when we obtain the 
radiative metric, we shall denote the associated radiative coordinates by (T,R). 
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where k a = (l,n) is the outgoing Minkowskian null vector. Given any n > 2, let us recursively 
assume that we have obtained all the previous radiative post-Minkowskian coefficients H"f 7 i.e. 
Vm < n — 1, and that all of them satisfy 

k^H^ = O (Jz) . (54) 

From this induction hypothesis one can prove that the nth post-Minkowskian source term A" 13 = 
A^(H 1 , . . . , ff„_i) is such that 

A n P = —Vn(u 7 n) + 0l^-y (55) 

To the leading order this term takes the classic form of the stress-energy tensor for a swarm 
of massless particles, with tr„ being related to the power in the waves. One can show that all 
the problems with the appearance of logarithms come from the retarded integral of the terms in 
Equation l|55|l that behave like 1/r 2 : See indeed the integration formula l|109|l . which behaves like 
In r/r at infinity. But now, thanks to the particular index structure of the term l|55|l . we can 
correct for the effect by adjusting the gauge at the nth post-Minkowskian order. We pose, as a 
gauge vector, 

— J dva n (v,n) , (56) 

where TV refers to the same finite part operation as in Equation l|39fl ■ This vector is such that the 
logarithms that will appear in the corresponding gauge terms cancel out the logarithms coming 
from the retarded integral of the source term Q55J1: see Ref. f° r the details. Hence, to the nth 
post-Minkowskian order, we define the radiative metric as 

Hf = Uf + + d a C n + - v af3 d^- (57) 

Here U^ 3 and V" 13 denote the quantities that are the analogues of u^j 3 and v" 13 , which were 
introduced into the harmonic-coordinates algorithm: See Equations (|39l I4UI HTl I42|) . In particular, 
these quantities are constructed in such a way that the sum Z7"^ + V"^ is divergence-free, so we 
see that the radiative metric does not obey the harmonic-gauge condition: 



C = Tvn 



ret 



d„H^ = n$ = ^ fdva n (v,n) 



(58) 



The far-zone expansion of the latter metric is of the type (|49f) . i.e. is free of any logarithms, and the 
retarded time in these coordinates tends asymptotically toward a null coordinate at infinity. The 
property of asymptotic simplicity, in the mathematical form given by Geroch and Horowitz |118) . is 
proved by introducing the conformal factor f2 = 1/r in radiative coordinates (see Ref. ^2]). Finally, 
it can be checked that the metric so constructed, which is a functional of the source multipole 
moments 1^ , . . . , (from the definition of the algorithm) , is as general as the general harmonic- 
coordinate metric of Theorem [21 since it merely differs from it by a coordinate transformation 
(t, x) — ► (T, X), where (t, x) are the harmonic coordinates and (T, X) the radiative ones, together 
with a re-definition of the multipole moments. 

4.4 The radiative multipole moments 

The leading-order term 1/R of the metric in radiative coordinates, neglecting 0(1/ R 2 ), yields 
the operational definition of two sets of STF radiative multipole moments, mass-type XJl(U) and 
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current-type Vl{U). By definition, we have 



AG + °° if 21 1 

ffJ T (t/, N) = -^V, jab (ti) Y, J5J |Wi- 3 U„n_ 2 (I/) - jjy— jJV oL _ 2e „, ( „V, )di _ 2 (i;) j 

+ (59) 

This multipole decomposition represents the generalization, up to any post-Newtonian order (wit- 
ness the factors of 1/c in front of each of the multipolar pieces) of the quadrupole-moment formalism 
reviewed in Equation @ ■ The corresponding total gravitational flux reads 

- G J (i + lKi + 2) TT (i) , mTT d)rmx 4 ^ + 2 ) 



;=2 



£ ^ - > ^^ (M).^!)!! ^^^ + ^(i - 1)(Z + l)!(2i + Dn V^C^V^Ct/) 



(60) 



Notice that the meaning of such formulas is rather empty, because we do not know yet how the 
radiative moments are given in terms of the actual source parameters. Only at the Newtonian 
level do we know this relation, which from the comparison with the quadrupole formalism of 
Equations (0 |21 reduces to 

V ij (U) = Q%\u) + o(±) ) (61) 



where Qy is the Newtonian quadrupole given by Equation J3J . Fortunately, we are not in such bad 
shape because we have learned from Theorem 0] the general method that permits us to compute 
the radiative multipole moments Ul, Vl in terms of the source moments II, Jl, . . . , Z^. Therefore, 
what is missing is the explicit dependence of the source multipole moments as functions of the actual 
parameters of some isolated source. We come to grips with this question in the next Section [S] 
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5 Exterior Field of a Post-Newtonian Source 



By Theorem [21 we control the most general class of solutions of the vacuum equations outside 
the source, in the form of non-linear functionals of the source multipole moments. For instance, 
these solutions include the Schwarzschild and Kerr solutions, as well as all their perturbations. 
By Theorem 0] we learned how to construct the radiative moments at infinity. We now want to 
understand how a specific choice of stress-energy tensor T a @ (i.e. a choice of some physical model 
describing the source) selects a particular physical exterior solution among our general class. 

5.1 The matching equation 

We shall provide the answer in the case of a post-Newtonian source for which the post-Newtonian 
parameter 1/c defined by Equation is small. The fundamental fact that permits the connection 
of the exterior field to the inner field of the source is the existence of a "matching" region, in which 
both the multipole and the post-Newtonian expansions are valid. This region is nothing but the 
exterior near zone, such that r > a (exterior) and r <C A (near zone). It always exists around 
post-Newtonian sources. 

Let us denote by A4(h) the multipole expansion of h (for simplicity, we suppress the space-time 
indices). By A4(h) we really mean the multipolar-post-Minkowskian exterior metric that we have 
constructed in Sections |21 and 

+oo 

M(h) = h cxt = Y, G n h n [l L , . . . , Z4 (62) 

71=1 

Of course, h agrees with its own multipole expansion in the exterior of the source, 

r>a =>■ M(h) = h. (63) 

By contrast, inside the source, h and A4(h) disagree with each other because h is a fully-fledged 
solution of the field equations with matter source, while M(h) is a vacuum solution becoming 
singular at r = 0. Now let us denote by h the post-Newtonian expansion of h. We have already 
anticipated the general structure of this expansion as given in Equation l|47[) . In the matching 
region, where both the multipolar and post-Newtonian expansions are valid, we write the numerical 
equality 

a<r<A =>■ M(h) = h. (64) 

This "numerical" equality is viewed here in a sense of formal expansions, as we do not control the 
convergence of the series. In fact, we should be aware that such an equality, though quite natural 
and even physically obvious, is probably not really justified within the approximation scheme 
(mathematically speaking), and we take it as part of our fundamental assumptions. 

We now transform Equation (|64|) into a matching equation, by replacing in the left-hand side 
A4(h) by its near-zone re-expansion J\A(h), and in the right-hand side h by its multipole expansion 
Ai(h). The structure of the near-zone expansion (r — > 0) of the exterior multipolar field has been 
found in Equation (|46J1 . We denote the corresponding infinite series Ai(h) with the same overbar 
as for the post-Newtonian expansion because it is really an expansion when r/c — > 0, equivalent 
to an expansion when c — > oo. Concerning the multipole expansion of the post-Newtonian metric, 
M.(h), we simply postulate its existence. Therefore, the matching equation is the statement that 

M(K)=M(R), (65) 

by which we really mean an infinite set of functional identities, valid V(x, t) £ x R, between 
the coefficients of the series in both sides of the equation. Note that such a meaning is somewhat 
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different from that of a numerical equality like Equation H64|) . which is valid only when x belongs 
to some limited spatial domain. The matching equation (|65|) tells us that the formal near-zone 
expansion of the multipole decomposition is identical, term by term, to the multipole expansion 
of the post-Newtonian solution. However, the former expansion is nothing but the formal far- 
zone expansion, when r — > oo, of each of the post-Newtonian coefficients. Most importantly, it 
is possible to write down, within the present formalism, the general structure of these identical 
expansions as a consequence of Theorem Equation l|46|l : 

A*W=$>Lr m (lnr)PF iimiP (t) = M(h), (66) 

where the functions Fi,,m,p — S n >i ^ n ^L,m,p,n- The latter expansion can be interpreted ei- 
ther as the singular re-expansion of the multipole decomposition when r — > (first equality in 
Equation Ij66|l ). or the singular re-expansion of the post-Newtonian series when r — > +oo (sec- 
ond equality). We recognize the beauty of singular perturbation theory where two asymptotic 
expansions, taken formally outside their respective domains of validity, are matched together. Of 
course, the method works because there exists, physically, an overlapping region in which the two 
approximation series are expected to be numerically close to the exact solution. 

5.2 General expression of the multipole expansion 

Theorem 5 Under the hypothesis of matching, Equation |6'5j) . the multipole expansion of the 
solution of the Einstein field equation outside a post-Newtonian source reads 

M(h^) = TV B=0 a; c l[¥ B M(A^)} ^^^tldJ-Hfit - r/c)\ , (67) 

c i=o ^ r J 

where the "multipole moments" are given by 

nf(u) = TV B =of d 3 x|xJ s a; L T Q/3 (x,u). (68) 

Here, t q/3 denotes the post- Newtonian expansion of the stress-energy pseudo-tensor defined by 
Equation \14\ l- 

Proof |15l I2U| : First notice where the physical restriction of considering a post-Newtonian source 
enters this theorem: the multipole moments l|68|l depend on the post-Newtonian expansion r"' 3 , 
rather than on r a P itself. Consider A"' 3 , namely the difference between ft"' 3 , which is a solution of 
the field equations everywhere inside and outside the source, and the first term in Equation l|67|) . 
namely the finite part of the retarded integral of the multipole expansion A4(A a @): 

A a/3 = h *f3 _ TV □- t i[^(A Q/3 )]. (69) 

From now on we shall generally abbreviate the symbols concerning the finite-part operation at 
B = by a mere TV . According to Equation (|20|) . h a P is given by the retarded integral of the 
pseudo-tensor r Q/3 . So, 

A «* = 1*1° - TV n£ [M(A«*)] . (70) 

In the second term the finite part plays a crucial role because the multipole expansion A4(A a/3 ) 
is singular at r = 0. By contrast, the first term in Equation JTUJl, as it stands, is well-defined 
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because we are considering only some smooth field distributions: T af} £ C°°(R 4 ). There is no need 
to include a finite part TV in the first term, but a contrario there is no harm to add one in front of 
it, because for convergent integrals the finite part simply gives back the value of the integral. The 
advantage of adding "artificially" the TV in the first term is that we can re-write Equation (|70() 
into the much more interesting form 

A<" = ^TV n-l [r^ - M(r^)] , (71) 

in which wc have also used the fact that M(A a P) = 16ttG/c 4 • M{r al3 ) because T a/3 has a compact 
support. The interesting point about Equation 171f) is that A a @ appears now to be the (finite part 
of a) retarded integral of a source with spatially compact support. This follows from the fact that 
the pseudo-tensor agrees numerically with its own multipole expansion when r > a (same equation 
as JHIIJ)). Therefore, Ai(A a/3 ) can be obtained from the known formula for the multipole expansion 
of the retarded solution of a wave equation with compact-support source. This formula, given in 
Appendix B of Ref. [2E] , yields the second term in Equation i|t)7ll , 



1=0 



but in which the moments do not yet match the result l|t)8[) : instead, 

Uf = TV J d 3 xx L [r a P - M{T al3 )] . (73) 

The reason is that we have not yet applied the assumption of a post-Newtonian source. Such 
sources are entirely covered by their own near zone (i.e. a <C A), and, in addition, the integral (|73l 
has a compact support limited to the domain of the source. In consequence, we can replace the 
integrand in Equation (|73|) by its post-Newtonian expansion, valid over all the near zone, i.e. 



Uf =TV I d 3 xx L 



^-M(t^) . (74) 



Strangely enough, we do not get the expected result because of the presence of the second term in 
Equation i|74|l . Actually, this term is a bit curious, because the object Ai(r al3 ) it contains is only 
known in the form of the formal series whose structure is given by the first equality in Equation 
(indeed r and h have the same type of structure). Happily (because we would not know what to do 
with this term in applications), we are now going to prove that the second term in Equation (|74(l 
is in fact identically zero. The proof is based on the properties of the analytic continuation as 
applied to the formal structure l|66|l of Ai(T a/3 ). Each term of this series yields a contribution 
to Equation H74J) that takes the form, after performing the angular integration, of the integral 
TVb=o Jq °° drr B+b (In r) p , and multiplied by some function of time. We want to prove that the 
radial integral J +O ° dr r B+b (In r) p is zero by analytic continuation (V-B £ C). First we can get rid 
of the logarithms by considering some repeated differentiations with respect to B; thus we need 
only to consider the simpler integral / +o ° drr B+h . We split the integral into a "near-zone" integral 

drr B+b and a "far-zone" one J^°° dr r B+b , where 1Z is some constant radius. When 3?(-B) is a 
large enough positive number, the value of the near-zone integral is lZ B+b+1 /(B + b+1), while when 
"St (B) is a large negative number, the far-zone integral reads the opposite, — lZ B+b+1 / (B + 6 + 1). 
Both obtained values represent the unique analytic continuations of the near-zone and far-zone 
integrals for any BeC except —6—1. The complete integral / +o ° drr B+b is equal to the sum of 
these analytic continuations, and is therefore identically zero (V-B £ C, including the value —6—1). 
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At last we have completed the proof of Theorem 

Hf =TV J d 3 xx L T aP . (75) 

The latter proof makes it clear how crucial the analytic-continuation finite part TV is, which 
we recall is the same as in our iteration of the exterior post-Minkowskian field (see Equation l|39[) 1. 
Without a finite part, the multipole moment (|75(l would be strongly divergent, because the pseudo- 
tensor t 01 ® has a non-compact support owing to the contribution of the gravitational field, and the 
multipolar factor xl behaves like r l when r — ► +00. In applications (Part [H] of this article) we 
must carefully follow the rules for handling the TV operator. 

The two terms in the right-hand side of Equation (|67|l depend separately on the length scale 
ro that we have introduced into the definition of the finite part, through the analytic-continuation 
factor r B = (r/r ) B (see Equation However, the sum of these two terms, i.e. the exterior 

multipolar field A4(h) itself, is independent of ro- To see this, the simplest way is to differentiate 
formally M{h) with respect to ro. The independence of the field upon ro is quite useful in ap- 
plications, since in general many intermediate calculations do depend on ro, and only in the final 
stage does the cancellation of the ro's occur. For instance, we shall see that the source quadrupolc 
moment depends on ro starting from the 3PN level but that this ro is compensated by another 
ro coming from the non-linear "tails of tails" at the 3PN order. 



5.3 Equivalence with the Will— Wiseman formalism 

Recently, Will and Wiseman |217| (see also Refs. |216l 1170) 1. extending previous work of Epstein 
and Wagoner |1()7] and Thorne |2()7| . have obtained a different-looking multipole decomposition, 
with different definitions for the multipole moments of a post-Newtonian source. They find, instead 
of our multipole decomposition given by Equation (|fJ7)l , 

M(h^)=U- c \[M{K a % - —Y l L± dL W f(i-r/c) . (76) 

l TC c ;=o t. ^r j 

There is no TV operation in the first term, but instead the retarded integral is truncated, as 
indicated by the subscript 1Z, to extend only in the "far zone": i.e. |x'| > 1Z in the notation of 
Equation l|21|l . where TZ is a constant radius enclosing the source {TZ > a). The near-zone part 
of the retarded integral is thereby removed, and there is no problem with the singularity of the 
multipole expansion Ai(A a @) at the origin. The multipole moments Wl arc then given, in contrast 
with our result I|tj8fl . by an integral extending over the "near zone" only: 

Wf(u)= f d^x L T a P{^u). (77) 

Since the integrand is compact-supported there is no problem with the bound at infinity and the 
integral is well-defined (no need of a TV). 

Let us show that the two different formalisms are equivalent. We compute the difference 
between our moment TLl, defined by Equation l|68[l . and the Will- Wiseman moment Wl, given by 
Equation (|77(l . For the comparison we split TLl into far-zone and near-zone integrals corresponding 
to the radius TZ. Since the finite part TV present in Ti^ deals only with the bound at infinity, it 
can be removed from the near-zone integral, which is then seen to be exactly equal to Wl- So the 
difference between the two moments is simply given by the far-zone integral: 

n a L p {u)-Wf{u)=Tv( d^XL^i^u). (78) 

J\x\>TZ 
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Next, we transform this expression. Successively we write r a ^ = M{r a ^) because we are outside 
the source, and M(j al3 ) = M.(r a P) from the matching equation |[BSJ. At this stage, we recall from 
our reasoning right after Equation (|74|) that the finite part of an integral over the whole space M 3 
of a quantity having the same structure as M(T a @) is identically zero by analytic continuation. 
The main trick of the proof is made possible by this fact, as it allows us to transform the far-zone 
integration |x| > 1Z in Equation (|78(l into a near-zone one |x| < 7?., at the price of changing the 
overall sign in front of the integral. So, 

nf(u) - Wf(u) = -TV f d 3 X x L M(T<*0)( X ,u). (79) 

J |x|<7Z 

Finally, it is straightforward to check that the right-hand side of this equation, when summed 
up over all multipolarities I, accounts exactly for the near-zone part that was removed from the 
retarded integral of A4(A a P) (first term in Equation ((7SJ), so that the "complete" retarded integral 
as given by the first term in our own definition l|67[) is exactly reconstituted. In conclusion, the 
formalism of Ref. |217| is equivalent to the one of Refs. |15II2(J) . 



5.4 The source multipole moments 

In principle the bridge between the exterior gravitational field generated by the post-Newtonian 
source and its inner field is provided by Theorem |3J however, we still have to make the con- 
nection with the explicit construction of the general multipolar and post-Minkowskian metric in 
Sections [3] and Namely, we must find the expressions of the six STF source multipole moments 
II, Jl, ■ ■ ■ , Zl parametrizing the linearized metric (|26l I27ll2"5|l at the basis of that construction 12 . 

To do this we first find the equivalent of the multipole expansion given in Theorem^] which was 
parametrized by non-trace-free multipole functions Tt^( u ), m terms of new multipole functions 
rf \u) that are STF in all their indices L. The result (which follows from Equation (B.14a) 
in |2EJ) is 

Miff*) = TV n-l[M(A<*?)} - — £ L± dL Tf{t - r/c) , (80) 

i=o " ^ ' 

where the STF multipole functions (witness the multipolar factor xl = STF[xl]) read 
Tf(u)=TV J d^xL J < 



dz5i{z)T an (x,u + z\x\/c). (81) 

Notice the presence of an extra integration variable z, ranging from —1 to 1. The z-integration 
involves the weighting function 13 

S^^^^il-z 2 ) 1 , (82) 
which is normalized in such a way that 

r 1 

dzSi(z) = l. (83) 



i-i 

12 Recall that in actual applications we need mostly the mass-type moment I/, and current-type one J^, because 
the other moments parametrize a linearized gauge transformation. 

13 This function approaches the Dirac delta-function (hence its name) in the limit of large multipolcs: 
limj^+oo <5;(z) = S(z). Indeed the source looks more and more like a point mass as we increase the multipolar 
order I. 
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The next step is to impose the harmonic-gauge conditions 1|12|) onto the multipole decomposi- 
tion JSUJ, and to decompose the multipole functions T^f {u) into STF irreducible pieces with 
respect to both L and their space-time indices a/3. This technical part of the calculation is identi- 
cal to the one of the STF irreducible multipole moments of linearized gravity [EH]. The formulas 
needed in this decomposition read 



ttW _ ( + ) T , (0) T +5^ (-^Tr I--. 

•> L — J-iL \ &ai<ii ± L — l>a ^ u i<%i ± L— 1>; 

(84) 

J* = ^U l]L + STF STF [e aill {+1) U a3L ^ + W^x-r 

L 13 

+S tll e ajll _ 1 (-VU aL -2 + S tll 5 j i [ _ 1 ( --^U L - 2 } + S lj V L , 

where the ten tensors i?L, ^ + 'Tl+i, • • ■ , '~ 2 ^l-2, Vl are STF, and are uniquely given in terms 
of the T\ 
theorem. 



of the Tl p, s by some inverse formulas. Finally, the latter decompositions lead to the following 



Theorem 6 The STF multipole moments I/, and Jl of a post-Newtonian source are given, formally 
up to any post- Newtonian order, by (I > 2) 



I £ («) = TV J *x £ dz{s^ - g^ ^ 1 ' 



2(2/ + 1) r A ^(2) 

c 4 (Z + l)(Z + 2)(2Z + 5) 



5 l+2 x ijL X\? }>(x,u + z|x|/c), (85) 



J L (u)=FV J d 3 x J dz e ab{il |<5;x L _i) a S & - ^ ^TjjT^ + ^ ^+i^L~i) ac S^ ) | (x, u + z|x|/c). 

These moments are the ones that are to be inserted into the linearized metric h°[ that represents 
the lowest approximation to the post-Minkowskian field h"^ t — 2 n >i G n h"@ defined in Section^ 

In these formulas the notation is as follows: Some convenient source densities are defined from the 
post-Newtonian expansion of the pseudo-tensor T a @ by 

-00 + -„: 

v _ r 01 (86) 



T ■ ■ — ¥ l - 



(where r M = Siff 1 ^). As indicated in Equations l|85|l these quantities are to be evaluated at the 
spatial point x and at time u + z|x|/c. 

For completeness, we give also the formulas for the four auxiliary source moments Wl, . . . , Z^, 
which parametrize the gauge vector <pf as defined in Equations l|28|l : 



2 *+! , . yd) 



2c 2 (Z + 1)(Z + 2)(2Z + 5) 



(87) 
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x L (u) = rpj**^ d*{ 2{l + r^rm + 5 ) Si+2i ^}' 



-i 



Yl(u) = TV Id x / dzl -diXL^u + 7-— -p-7-— — ■5; + ii li E 



3(2Z + 1) , . 



(/ + l)(2I + 3) 
2(2/ + 1) ( a)\ 

r d i+ 2^yxS- 5 - h (89) 



Zi(u) =TV I d 3 x / dze ab(4l <j -— — — — ^— <5; + ii L _i )6c £ ac }>. (90) 



c 2 (/ + l)(/ + 2)(2; + 5) 



2Z + 1 



1 



(Z + 2)(2Z + 3) 



As discussed in Section^] one can always find two intermediate "packages" of multipole moments, 
Ml and S l , which are some non-linear functionals of the source moments l|85ll and Equations 1)871 
IHHIIHSlinnij and such that the exterior field depends only on them, modulo a change of coordinates 
(see, e.g., Equation l(55|l below). 

In fact, all these source moments make sense only in the form of a post-Newtonian expansion, 
so in practice we need to know how to expand all the z- integrals as series when c — > +00. Here is 
the appropriate formula: 

\\ dZ MZ)T(X ' U + MIC) = I 2*^ + 2^1)!! (rif^ «)• < 91 ) 

Since the right-hand side involves only even powers of 1/c, the same result holds equally well for 
the "advanced" variable u + z\x\/c or the "retarded" one u — z\x\/c. Of course, in the Newtonian 
limit, the moments If, and Jl (and also Mi, Sl) reduce to the standard expressions. For instance, 
we have 

Il(u) = Q l (u) + o(^), (92) 

where is the Newtonian mass- type multipole moment (see Equation ©)■ (The moments 
Wi, . . . , lih have also a Newtonian limit, but it is not particularly illuminating.) 

Needless to say, the formalism becomes prohibitively difficult to apply at very high post- 
Newtonian approximations. Some post-Newtonian order being given, we must first compute the 
relevant relativistic corrections to the pseudo stress-energy-tensor r Q ^ (this necessitates solving 
the field equations inside the matter, see Section 15. 5f) before inserting them into the source mo- 
ments P51 1551 IS2l I5H l9"Tl 1571 1551 1551 150)l . The formula (JSJ is used to express all the terms up to 
that post-Newtonian order by means of more tractable integrals extending over R 3 . Given a spe- 
cific model for the matter source we then have to find a way to compute all these spatial integrals 
(we do it in Section fTTTl in the case of point-mass binaries). Next, we must substitute the source 
multipole moments into the linearized metric 1)261 1271 1281) , and iterate them until all the necessary 
multipole interactions taking place in the radiative moments Ul and Vl are under control. In fact, 
we shall work out these multipole interactions for general sources in the next section up to the 3PN 
order. Only at this point does one have the physical radiation field at infinity, from which we can 
build the templates for the detection and analysis of gravitational waves. We advocate here that 
the complexity of the formalism reflects simply the complexity of the Einstein field equations. It is 
probably impossible to devise a different formalism, valid for general sources devoid of symmetries, 
that would be substantially simpler. 
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5.5 Post-Newtonian field in the near zone 

Theorem \R\ solves in principle the question of the generation of gravitational waves by extended 
post-Newtonian sources. However, note that this result has to be completed by the definition 
of an explicit algorithm for the post-Newtonian iteration, analogous to the post-Minkowskian al- 
gorithm we defined in Section 01 so that the source multipole moments, which contain the full 
post-Newtonian expansion of the pseudo-tensor r"' 3 , can be completely specified. Such a sys- 
tematic post-Newtonian iteration scheme, valid (formally) to any post-Newtonian order, has been 
implemented [T%2"1 EH] using matched asymptotic expansions. The solution of this problem yields, 
in particular, some general expression, valid up to any order, of the terms associated with the 
gravitational radiation reaction force inside the post-Newtonian source . 

Before proceeding, let us recall that the "standard" post-Newtonian approximation, as it was 
used until, say, the early 1980's (see for instance Refs. |2l 11391 IT4"U1 1169) 1. is plagued with some 
apparently inherent difficulties, which crop up at some high post-Newtonian order. The first 
problem is that in higher approximations some divergent Poisson-type integrals appear. Indeed 
the post-Newtonian expansion replaces the resolution of a hyperbolic-like d'Alembertian equation 
by a perturbatively equivalent hierarchy of elliptic-like Poisson equations. Rapidly it is found 
during the post-Newtonian iteration that the right-hand side of the Poisson equations acquires a 
non-compact support (it is distributed over all space), and that as a result the standard Poisson 
integral diverges at the bound of the integral at spatial infinity, i.e. r = |x| — > +oo, with t — const. 

The second problem is related with the a priori limitation of the approximation to the near 
zone, which is the region surrounding the source of small extent with respect to the wavelength 
of the emitted radiation: r <C A. The post-Newtonian expansion assumes from the start that all 
retardations r/c are small, so it can rightly be viewed as a formal near-zone expansion, when r — * 0. 
In particular, the fact which makes the Poisson integrals to become typically divergent, namely 
that the coefficients of the post-Newtonian series blow up at "spatial infinity", when r — > +oo, has 
nothing to do with the actual behaviour of the field at infinity. However, the serious consequence 
is that it is not possible, a priori, to implement within the post-Newtonian iteration the physical 
information that the matter system is isolated from the rest of the universe. Most importantly, 
the no-incoming radiation condition, imposed at past null infinity, cannot be taken into account, 
a priori, into the scheme. In a sense the post- Newtonian approximation is not "self-supporting" , 
because it necessitates some information taken from outside its own domain of validity. 

Here we present, following Refs. [1821 141) . a solution of both problems, in the form of a general 
expression for the near-zone gravitational field, developed to any post-Newtonian order, which 
has been determined from implementing the matching equation This solution is free of 

the divergences of Poisson-type integrals we mentionned above, and it incorporates the effects of 
gravitational radiation reaction appropriate to an isolated system. 

Theorem 7 The expression of the post-Newtonian field in the near zone of a post-Newtonian 
source, satisfying correct boundary conditions at infinity (no incoming radiation), reads 



—a(3 16-kG 

tl = ; 



a \nf{t-r/c)-nf{t + r/c) 



Z! 2r 



(93) 



The first term represents a particular solution of the hierarchy of post-Newtonian equations, while 
the second one is a homogeneous multipolar solution of the wave equation, of the "anti-symmetric" 
type that is regular at the origin r — located in the source. 



14 An alternative approach to the problem of radiation reaction, besides the matching procedure, is to work only 
within a post-Minkowskian iteration scheme (which does not expand the retardations): see, e.g., Ref. 1691 . 
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More precisely, the flat retarded d'Alembertian operator in Equation ((93(1 is given by the standard 
expression 1(21(1 but with all retardations expanded (r/c — > 0), and with the finite part J-V proce- 
dure involved for dealing with the bound at infinity of the Poisson-type integrals (so that all the 
integrals are well-defined at any order of approximation), 

rpu£ [t**] = -i- £ tf- ( AV ^ / d 3 x , |x _ xT ,-i^ (x ^). (94) 

The existence of the solution 194(1 shows that the problem of divergences of the post-Newtonian 
expansion is simply due to the fact that the standard Poisson integral does not constitute the 
correct solution of the Poisson equation in the context of post-Newtonian expansions. So the 
problem is purely of a technical nature, and is solved once we succeed in finding the appropriate 
solution to the Poisson equation. 

Theorem is furthermore to be completed by the information concerning the multipolar func- 
tions 1Z™ (it) parametrizing the anti-symmetric homogeneus solution, the second term of Equa- 
tion 1(93(1 . Note that this homogeneous solution represents the unique one for which the matching 
equation ((d3)| is satisfied. The result is 

1 f r+oo 

Kf(u) = -^J r Vjd 3 x.'x' L J d zll {z)M{T al3 ){^,u-z\^\lc), (95) 

where M. (r a P ) denotes the multipole expansion of the pseudo-tensor (in the sense of Equation l(62(l ) , 
and where we denote 7/(2) = — 2<5;(z), with 6i(z) being given by Equation ((82(1 15 . 

Importantly, we find that the post-Newtonian expansion h given by Theorem0is a functional 
not only of the related expansion of the pseudo-tensor, r Q/3 , but also, by Equation lj9*5)l. of its 
multipole expansion M(r°' ? ), which is valid in the exterior of the source, and in particular in 
the asymptotic regions far from the source. This can be understood by the fact that the post- 
Newtonian solution ((93(1 depends on the boundary conditions imposed at infinity, that describe a 
matter system isolated from the rest of the universe. 

Equation ((93|) is interesting for providing a practical recipe for performing the post-Newtonian 
iteration ad infinitum. Moreover, it gives some insights on the structure of radiation reaction terms. 
Recall that the anti-symmetric waves, regular in the source, are associated with radiation reaction 
effects. More precisely, it has been shown 182 that the specific anti-symmetric wave given by 
the second term of Equation 1(93(1 is linked with some non-linear contribution due to gravitational 
wave tails in the radiation reaction force. Such a contribution constitutes a generalization of 
the tail-transported radiation reaction term at the 4PN order, i.e. 1.5PN order relative to the 
dominant radiation reaction order, as determined in Ref. |27) . This term is in fact required by 
energy conservation and the presence of tails in the wave zone (see, e.g., Equation (|§7)l below). 
Hence, the second term of Equation ((93(1 is dominantly of order 4PN and can be neglected in 
computations of the radiation reaction up to 3.5PN order (as in Ref. (161( 1. The usual radiation 
reaction terms, up to 3.5PN order, which are linear in the source multipole moments (for instance 
the usual radiation reaction term at 2.5PN order), are contained in the first term of Equation 1(93(1 . 
and are given by the terms with odd powers of 1 / c in the post-Newtonian expansion ((94|) . It can 
be shown that such terms take also the form of some anti-symmetric multipolar wave, which 
turn out to be parametrized by the same moments as in the exterior field, namely the moments 
which are the STF analogues of Equations ((68|) . 



15 Notice that the normalization J^°° dzji(z) = 1 holds as a consequence of the corresponding normalization 1831 
for <5;(z), together with the fact that dz 7; (2) = by analytic continuation in the variable I £ C. 
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6 Non-linear Multipole Interactions 



We shall now show that the radiative mass-type quadrupole moment U y - includes a quadratic 
tail at the relative 1.5PN order (or 1/c 3 ), corresponding to the interaction of the mass M of the 
source and its quadrupole moment Iy. This is due to the back-scattering of quadrupolar waves 
off the Schwarzschild curvature generated by M. Next, Uy includes a so-called non-linear memory 
integral at the 2.5PN order, due to the quadrupolar radiation of the stress-energy distribution of 
linear quadrupole waves themselves, i.e. of multipole interactions I y - x Iki- Finally, we have also a 
cubic tail, or "tail of tail" , arising at the 3PN order, and associated with the multipole interaction 
r2 x Ijj . The result for U y - is better expressed in terms of the intermediate quadrupole moment 



M 2 



My already discussed in Section l4~2l This moment reads |16) 



I 

C J 



(96) 



where W means Wl as given by Equation (|87() in the case / = (of course, in Equation (|96|l we 
need only the Newtonian value of W). The difference between the two moments and ly is 
a small 2.5PN quantity. Henceforth, we shall express many of the results in terms of the mass 
moments and the corresponding current ones S^. The complete formula for the radiative 
quadrupole, valid through the 3PN order, reads [21 EH 
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(97) 



The retarded time in radiative coordinates is denoted U = T — R/c. The constant tq is the one 
that enters our definition of the finite-part operation TV (see Equation (|36|l '). The "Newtonian" 
term in Equation (|9"?|l contains the Newtonian quadrupole moment Qy (see Equation (|9"2"|0 . The 
dominant radiation tail at the 1.5PN order was computed within the present formalism in Ref. |29j . 
The 2.5PN non-linear memory integral - the first term inside the coefficient of G/c 5 - has been 
obtained using both post-Newtonian methods ^3 12191 12101 l2"9"l l2"T] and rigorous studies of the 
field at future null infinity [7J. The other multipole interactions at the 2.5PN order can be found 
in Ref. |3T]. Finally the "tail of tail" integral appearing at the 3PN order has been derived in 
this formalism in Ref. 19. ■ Be careful to note that the latter post-Newtonian orders correspond 
to "relative" orders when counted in the local radiation-reaction force, present in the equations 
of motion: For instance, the 1.5PN tail integral in Equation (|97|) is due to a 4PN radiative effect 
in the equations of motion [27]; similarly, the 3PN tail-of-tail integral is (presumably) associated 
with some radiation-reaction terms occuring at the 5.5PN order. 

Notice that all the radiative multipole moments, for any I, get some tail- induced contributions. 
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They are computed at the 1.5PN level in Appendix C of Ref. ^H]- We find 



V L (U)=M ( 1\U) + 
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(98) 



where the constants and 717 are given by 
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Recall that the retarded time U in radiative coordinates is given by 

2GM 



U = t 



''() 



In - + O (G 2 ) 



(100) 



where (t,r) are harmonic coordinates; recall the gauge vector £° in Equation l|51|l . Inserting U as 
given by Equation l|100fl into Equations (|98|) we obtain the radiative moments expressed in terms 
of source-rooted coordinates (i, r), e.g., 



U i = M«( i -r/c) + ^ 



dr M ( l +2) (t-T -r/c) 



In 



o 



(101) 



This expression no longer depends on the constant r$ (i.e. the vq gets replaced by r) 16 . If 
we now change the harmonic coordinates (t, r) to some new ones, such as, for instance, some 
"Schwarzschild-like" coordinates (t',r') such that t' — t and r' = r + GM/c 2 , we get 



U, 



-M^it' -r'/c) + 



2GM 



rfrM ( / +2) ^' 
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(102) 



where = ni + 1/2. Therefore the constant ki (and 7r; as well) depends on the choice of source- 
rooted coordinates (t, r): For instance, we have K2 = 11/12 in harmonic coordinates (see Equa- 
tion JH3), but k' 2 — 17/12 in Schwarzschild coordinates |50) . 

The tail integrals in Equations (|97l I98|) involve all the instants from —00 in the past up to the 
current time U. However, strictly speaking, the integrals must not extend up to minus infinity in 
the past, because we have assumed from the start that the metric is stationary before the date 
— T; see Equation i|19fl . The range of integration of the tails is therefore limited a priori to the 
time interval [-T, U]. But now, once we have derived the tail integrals, thanks in part to the 
technical assumption of stationarity in the past, we can argue that the results are in fact valid in 
more general situations for which the field has never been stationary. We have in mind the case 
of two bodies moving initially on some unbound (hyperbolic-like) orbit, and which capture each 
other, because of the loss of energy by gravitational radiation, to form a bound system at our 
current epoch. In this situation we can check, using a simple Newtonian model for the behaviour 

16 At the 3PN order (taking into account the tails of tails), we find that ro does not completely cancel out after 
the replacement of U by the right-hand side of Equation I1U0I . The reason is that the moment Mj, also depends 
on ro at the 3PN order. Considering also the latter dependence we can check that the 3PN radiative moment 
is actually free of the unphysical constant ro . 
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of the quadrupole moment My(f7 — r) when r — > +00, that the tail integrals, when assumed to 
extend over the whole time interval [—00, U], remain perfectly well-defined (i.e. convergent) at the 
integration bound r = +00. We regard this fact as a solid a posteriori justification (though not a 
proof) of our a priori too restrictive assumption of stationarity in the past. This assumption does 
not seem to yield any physical restriction on the applicability of the final formulas. 

To obtain the result (|97JI . we must implement in details the post-Minkows-kian algorithm 
presented in Section f4.il Let us outline here this computation, limiting ourselves to the interaction 
between one or two masses M = Madm = I and the time- varying quadrupole moment M (,(u) (that 
is related to the source quadrupole l a b(u) by Equation 196fl ). For these moments the linearized 
metric (J2E1 EH EHJ reads 

hf = h$ )+ h$ abV (103) 

where the monopole part is nothing but the linearized piece of the Schwarzschild metric in harmonic 
coordinates, 

h °(M) = -4r- x M, 

h °(M) = 0, (104) 



"(M) 
(M) 



and the quadrupole part is 



h°° = -2d ab [r^MM] 



(M 

hOi _ on 



•- x M<£>(u) , (105) 



(We pose c = 1 until the end of this section.) Consider next the quadratically non- linear metric 
h,2 generated by these moments. Evidently it involves a term proportional to M 2 , the mixed term 
corresponding to the interaction M x M a {,, and the self- interaction term of M a h. Say, 

a 2 (M 2 ) (MM ab ) ' (M ab M cd )' \ wo ) 

The first term represents the quadratic piece of the Schwarzschild metric, 

h °(M*) = 0, (107) 
^ J M2) = -n l3 r- 2 M 2 . 

The second term in Equation (|106|) represents the dominant non-static multipole interaction, that 
is between the mass and the quadrupole moment, and that we now compute . We apply Equa- 
tions (|39l 1401 mi 1421 143fl in Section 0] First we obtain the source for this term, viz. 

^(mm„i) = N af3 [h {M) , h {Mab) ] + N a P[h (Mab) , ( 108 ) 

where N a @(h } h) denotes the quadratic-order part of the gravitational source, as defined by Equa- 
tion H16JI . To integrate this term we need some explicit formulas for the retarded integral of an 
extended (non-compact-support) source having some definite multipolarity I. A thorough account 
of the technical formulas necessary for handling the quadratic and cubic interactions is given in 

17 The computation of the third term in Equation 11061 . which corresponds to the interaction between two 
quadrupoles, M ai) X M C[J , can be found in Ref. 1211 . 
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the appendices of Rcfs. and ^Hj- For the present computation the crucial formula corresponds 
to a source term behaving like 1/r 2 : 
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(109) 



where Qi is the Legendre function of the second kind 18 . With the help of this and other formulas 
we obtain the object given by Equation (|39() . Next we compute the divergence — d^u^, 
and obtain the supplementary term by applying Equations (|42|) . Actually, we find for this 
particular interaction w% — and thus also = 0- Following Equation (|43(l . the result is the 
sum of and , and we get 
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The metric is composed of two types of terms: "instantaneous" ones depending on the values of 
the quadrupole moment at the retarded time u = t — r, and "non-local" or tail integrals, depending 
on all previous instants t — rx < u. 

Let us investigate now the cubic interaction between two mass monopoles M with the quadrupole 
M Q f,. Obviously, the source term corresponding to this interaction reads 

^(M 2 M a |,) = Na>3 [h(M),h( MMab )} + iV a/? [fe( MMob ), fo(M)] + N al3 [h( M 2^ h(M ab )] + -W a/3 [/l(M al) ), /i(M 2 )] 



s The function Q { is given in terms of the Legendre polynomial P ; by 



1 r 1 dzP t (z) 1 



1 X — 2 
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In the complex plane there is a branch cut from — oo to 1. The first equality is known as the Neumann formula for 
the Legendre function. 
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■ M a/3 [h(M), h(M), h(u ab )] + M a(3 [h( M ), h( Mab ), ^(m)] + M 



a/3 



[^(M q6 ), h(M), /l(M)] 



(111) 



(see Equation l|33[l ). Notably, the TV-terms in Equation i|lll[l involve the interaction between 
a linearized metric, /i(M) or h(M ab ), and a quadratic one, /i(M 2 ) or /i(MM a6 )- So, included into 
these terms are the tails present in the quadratic metric /i(MM ab ) computed previously with the 
result IjllQp . These tails will produce in turn some "tails of tails" in the cubic metric ft-(M 2 M ab )- 
The rather involved computation will not be detailed here (see Ref. JH]). Let us just mention the 
most difficult of the needed integration formulas 19 : 
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where is the time anti-derivative of F. With this formula and others given in Ref. ^H] we 

are able to obtain the closed algebraic form of the metric ft^ 2M ^, at the leading order in the 
distance to the source. The net result is 
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where all the moments M Q b are evaluated at the instant t - 



r (recall that c = 1). Notice that 



some of the logarithms in Equations l|113fl contain the ratio r/r while others involve t/tq. The 



19 Equation 11121 has been obtained using a not so well known mathematical relation between the Legendre 
functions and polynomials: 

- / = Ql(x)Pl{y) 

2 7-i - 2)2 - (z 2 - l)( y 2 - 1) 

(where 1 < y < a; is assumed). See Appendix A in Ref. [191 for the proof. This relation constitutes a generalization 
of the Neumann formula (see footnote after Equation 11091 V 
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indicated remainders 0(l/r 2 ) contain some logarithms of r; in fact they should be more accurately 
written as o(r £_2 ) for some e<l. 

The presence of logarithms of r in Equations l|113[l is an artifact of the harmonic coordinates x a , 
and we need to gauge them away by introducing the radiative coordinates X a at future null infinity 
(see Theorem 0}. As it turns out, it is sufficient for the present calculation to take into account 
the "linearized" logarithmic deviation of the light cones in harmonic coordinates so that X a = 
x a +G£i +0(G 2 ), where £f is the gauge vector defined by Equation fBTJl (see also Equation l|10()jl ). 
With this coordinate change one removes all the logarithms of r in Equations (|113l) . Hence, we 
obtain the radiative metric 
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where the moments are evaluated at time U— r = T—R—t. It is trivial to compute the contribution 
of the radiative moments XJl(U) and Vl(U) corresponding to that metric. We find the "tail of 
tail" term reported in Equation H97fl . 



39 



7 The Third Post-Newtonian Metric 



The detailed calculations that are called for in applications necessitate having at one's disposal 
some explicit expressions of the metric coefficients g a p, in harmonic coordinates, at the highest 
possible post-Newtonian order. The 3PN metric that we present below 20 is expressed by means of 
some particular retarded-type potentials, V, Vi, Wij, etc., whose main advantages are to somewhat 
minimize the number of terms, so that even at the 3PN order the metric is still tractable, and to 
delineate the different problems associated with the computation of different categories of terms. 
Of course, these potentials have no physical significance by themselves. The basic idea in our 
post-Newtonian iteration is to use whenever possible a "direct" integration, with the help of some 



formulas like U'^d^Vd^V 
Ref. [SHI) reads 
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All the potentials are generated by the matter stress-energy tensor T a P through the definitions 
(analogous to Equations 1 (85) 1) 
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V and Vi represent some retarded versions of the Newtonian and gravitomagnetic potentials, 

V = D-\ [-4nGa] , 



Vi = □: 



(117) 



From the 2PN order we have the potentials 
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'Actually, such a metric is valid up to 3.5PN order. 
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Some parts of these potentials are directly generated by compact-support matter terms, while other 
parts are made of non-compact-support products of V-type potentials. There exists also a very 
important cubically non-linear term generated by the coupling between Wi 3 and V, the second 
term in the X-potential. At the 3PN level we have the most complicated of these potentials, 
namely 
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which involve many types of compact-support contributions, as well as quadratic-order and cubic- 
order parts; but, surprisingly, there are no quartically non-linear terms 21 . 

The above potentials are not independent. They are linked together by some differential iden- 
tities issued from the harmonic gauge conditions, which are equivalent, via the Bianchi identities, 
to the equations of motion of the matter fields (see Equation O). These identities read 
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It is important to remark that the above 3PN metric represents the inner post-Newtonian 
field of an isolated system, because it contains, to this order, the correct radiation-reaction terms 
corresponding to outgoing radiation. These terms come from the expansions of the retardations in 
the retarded- type potentials (fTTTl fTT^l fTTU)l . 



21 It has been possible to "integrate directly" all the quartic contributions in the 3PN metric. See the terms 
composed of V 4 and VX in the first of Equations 11151 . 
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Part II 

Compact Binary Systems 

The problem of the motion and gravitational radiation of compact objects in post-Newtonian ap- 
proximations of general relativity is of crucial importance, for at least three reasons. First, the 
motion of N objects at the 1PN level (1/c 2 ), according to the Einstein- Infold- Hoffmann equa- 
tions |1U6) . is routinely taken into account to describe the Solar System dynamics (see Ref. |16l)| ). 
Second, the gravitational radiation-reaction force, which appears in the equations of motion at the 
2.5PN order, has been experimentally verified, by the observation of the secular acceleration of the 
orbital motion of the binary pulsar PSR 1913+16 |2fl5l I50SI 1201) . 

Last but not least, the forthcoming detection and analysis of gravitational waves emitted by 
inspiralling compact binaries - two neutron stars or black holes driven into coalescence by emission 
of gravitational radiation - will necessitate the prior knowledge of the equations of motion and 
radiation field up to high post- Newtonian order. As discussed in the introduction in Section ^ 
(see around Equations @ El EJO) the appropriate theoretical description of inspiralling compact 
binaries is by two structureless point-particles, characterized solely by their masses mi and m-i 
(and possibly their spins), and moving on a quasi-circular orbit. Strategies to detect and analyze 
the very weak signals from compact binary inspiral involve matched filtering of a set of accurate 
theoretical template waveforms against the output of the detectors. Several analyses [TTJ EHJ lllll 
EH 12001 HHOI HS1 H3HI E2 E31 EHI EH EU lU El have shown that, in order to get sufficiently accurate 
theoretical templates, one must include post-Newtonian effects up to the 3PN level at least. 

To date, the templates have been completed through 3.5PN order for the phase evolution |351 
141)1 131| . and 2.5PN order for the amplitude corrections 001 ID Spin effects are known for the 
dominant relativistic spin-orbit coupling term at 1.5PN order and the spin-spin coupling term 
at 2PN order |143l El 11411 I7U) . and also for the next-to- leading spin-orbit coupling at 2.5PN 

order nnoi Enn inni En • 
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8 Regular izat ion of the Field of Point Particles 



Our aim is to compute the metric (and its gradient needed in the equations of motion) at the 
3PN order for a system of two point-like particles. A priori one is not allowed to use directly 
the metric expressions (|115|l . as they have been derived under the assumption of a continuous 
(smooth) matter distribution. Applying them to a system of point particles, we find that most of 
the integrals become divergent at the location of the particles, i.e. when x — > yx(t) or y2(i), where 
yi(t) and y2(i) denote the two trajectories. Consequently, we must supplement the calculation by 
a prescription for how to remove the "infinite part" of these integrals. At this stage different choices 
for a "self-field" regularization (which will take care of the infinite self-field of point particles) are 
possible. In this section we review: 

1. Hadamard's self-field regularization, which has proved to be very convenient for doing prac- 
tical computations (in particular, by computer), but suffers from the important drawback of 
yielding some ambiguity parameters, which cannot be determined within this regularization, 
at the 3PN order; 

2. Dimensional self-field regularization, an extremely powerful regularization which is free of any 
ambiguities (at least up to the 3PN level), and permits therefore to uniquely fix the values 
of the ambiguity parameters coming from Hadamard's regularization. However, dimensional 
regularization has not yet been implemented to the present problem in the general case (i.e. 
for an arbitrary space dimension d £ C) . 

The why and how the final results are unique and independent of the employed self-field regular- 
ization (in agreement with the physical expectation) stems from the effacing principle of general 
relativity |81| - namely that the internal structure of the compact bodies makes a contribution 
only at the formal 5PN approximation. However, we shall review several alternative computations, 
independent of the self-field regularization, which confirm the end results. 



8.1 Hadamard self-field regularization 

In most practical computations we employ the Hadamard regularization |125l 1196] (see Ref. |197| 
for an entry to the mathematical literature). Let us present here an account of this regularization, 
as well as a theory of generalized functions (or pseudo- functions) associated with it, following the 
investigations detailed in Refs. [3T>1 159) . 

Consider the class F of functions -F(x) which are smooth (C°°) on M 3 except for the two points 
yi and y2, around which they admit a power- like singular expansion of the type 22 

Vn S N, F(x) = r i {M + o(r?), (121) 

and similarly for the other point 2. Here r% = |x — yi| — > 0, and the coefficients i/ Q of the 
various powers of n depend on the unit direction ni = (x — yi)/r± of approach to the singular 
point. The powers a of n are real, range in discrete steps (i.e. a € (<2i)ieF})j and are bounded from 
below (do < a). The coefficients \f a (and 2 fa) for which a < can be referred to as the singular 
coefficients of F. If F and G belong to T so does the ordinary product FG, as well as the ordinary 
gradient diF. We define the Hadamard partie finie of F at the location of the point 1 where it is 
singular as 




(122) 



22 The function -F(x) depends also on time t, through for instance its dependence on the velocities vi(t) and V2(i), 
but the (coordinate) t time is purely "spectator" in the regularization process, and thus will not be indicated. 
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where dfli = dO(ni) denotes the solid angle element centered on yi and of direction ni. Notice 
that because of the angular integration in Equation (|122|) . the Hadamard partie finie is "non- 
distributive" in the sense that 

{FG)x ^ (F)i(G)i in general. (123) 

The non-distributivity of Hadamard's partie finie is the main source of the appearance of ambiguity 
parameters at the 3PN order, as discussed in Section 

The second notion of Hadamard partie finie (Pf) concerns that of the integral J <i 3 xF, which 
is generically divergent at the location of the two singular points yi and y2 (we assume that the 
integral converges at infinity). It is defined by 

Pf slS2 [d 3 xF= lim i f d 3 xF + 4n V ( —\ + 4tt In (— Vr 3 F) , + 1 <-► 2 I . 

(124) 

The first term integrates over a domain S(s) defined as R 3 from which the two spherical balls 
r\ < s and r 2 < s of radius s and centered on the two singularities, denoted 0(yi, s) and B(ya, s), 
are excised: <S(s) = R 3 \ S(yi, s) U 0(y2, s). The other terms, where the value of a function at 
point 1 takes the meaning (|122(l are such that they cancel out the divergent part of the first term 
in the limit where s — > (the symbol 1 <-> 2 means the same terms but corresponding to the other 
point 2). The Hadamard partie-finie integral depends on two strictly positive constants si and S2, 
associated with the logarithms present in Equation l|124fl . These constants will ultimately yield 
some gauge- type constants, denoted by r[ and r 2 , in the 3PN equations of motion and radiation 
field. See Ref. j3E] for alternative expressions of the partie-finie integral. 

We now come to a specific variant of Hadamard's regularization called the extended Hadamard 
regularization and defined in Refs. |36l I3H| . The basic idea is to associate to any Fe fa pseudo- 
function, called the partie finie pseudo-function Pf F, namely a linear form acting on functions G 
of F, and which is defined by the duality bracket 

VG g F, (Pf F, G) = Pf J d 3 xFG. (125) 

When restricted to the set V of smooth functions (i.e. G°°(IR 4 )) with compact support (obviously 
we have T> C F), the pseudo- function Pf F is a distribution in the sense of Schwartz |196| . The 
product of pseudo-functions coincides, by definition, with the ordinary pointwise product, namely 
PfF. Pf G = Pf(FG). In practical computations, we use an interesting pseudo-function, con- 
structed on the basis of the Riesz delta function |187j . which plays a role analogous to the Dirac 
measure in distribution theory, <$i(x) = 5(x — yi). This is the so-called delta-pseudo-function Pf Si 
defined by 

WFeF, (Pf<*i,F) =Pf Jd 3 x5 1 F= (F)i, (126) 

where (F)\ is the partie finie of F as given by Equation Ijl22|l . From the product of Pf 6\ with 
any Pf F we obtain the new pseudo- function Pf(F<5i), that is such that 

VG G F, (Pf(i^i), G) - (FG)i. (127) 

As a general rule, we are not allowed, in consequence of the "non-distributivity" of the Hadamard 
partie finie, Equation (|123|l . to replace F within the pseudo- function ¥i(Fb~i) by its regularized 
value: Pi(F6i) ^ {F)i Pf Si in general. It should be noticed that the object Pi(FSi) has no 
equivalent in distribution theory. 

Next, we treat the spatial derivative of a pseudo-function of the type PfF, namely 9j(Pf F). 
Essentially, we require (in Ref. [3H]) that the so-called rule of integration by parts holds. By 
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this we mean that we are allowed to freely operate by parts any duality bracket, with the all- 
integrated ("surface") terms always zero, as in the case of non-singular functions. This requirement 
is motivated by our will that a computation involving singular functions be as much as possible 
the same as if we were dealing with regular functions. Thus, by definition, 

VF,GeJ, (di{PtF),G) = -(di(FfG),F). (128) 

Furthermore, we assume that when all the singular coefficients of F vanish, the derivative of Pf F 
reduces to the ordinary derivative, i.e. di(Pf F) = Pf (diF). Then it is trivial to check that the 
rule H128JI contains as a particular case the standard definition of the distributional derivative [196) . 
Notably, we see that the integral of a gradient is always zero: (<9i(PfF),l) = 0. This should 
certainly be the case if we want to compute a quantity (e.g., a Hamiltonian density) which is 
defined only modulo a total divergence. We pose 

di(Pf F) = Pf (diF) + Di[F], (129) 

where Pf (diF) represents the "ordinary" derivative and Dj[F] the distributional term. The follow- 
ing solution of the basic relation l|128|) was obtained in Rcf. 36 : 



Di[F] =4tt Pf 



2 i ' rf i 

fe>0 1 



5i + 1 <-> 2, (130) 



where for simplicity we assume that the powers a in the expansion l|121fl of F are relative integers. 
The distributional term (|130l) is of the form Pf(G5i) (plus 1 <-> 2). It is generated solely by the 
singular coefficients of F (the sum over k in Equation (|130|) is always finite since there is a maximal 
order ao of divergency in Equation (|121|0 . The formula for the distributional term associated with 
the Ith distributional derivative, i.e. P>l[F] = 8l Pf F — Pf <9lF, where L = iii 2 . . . ii, reads 

i 

D L [F] » ... : IX /A f . (131) 

k=l 

We refer to Theorem 4 in Ref. j^H] for the definition of another derivative operator, representing 
the most general derivative satisfying the same properties as the one defined by Equation (|130|l . 
and, in addition, the commutation of successive derivatives (or Schwarz lemma) 23 . 

The distributional derivative (|129l 11301 1131|) does not satisfy the Leibniz rule for the derivation 
of a product, in accordance with a general result of Schwartz [1 951 . Rather, the investigation |36| 
suggests that, in order to construct a consistent theory (using the "ordinary" product for pseudo- 
functions), the Leibniz rule should be weakened, and replaced by the rule of integration by part, 
Equation (|128l) . which is in fact nothing but an "integrated" version of the Leibniz rule. However, 
the loss of the Leibniz rule stricto sensu constitutes one of the reasons for the appearance of the 
ambiguity parameters at 3PN order. 

The Hadamard regularization (F)i is defined by Equation (|122|) in a preferred spatial hyper- 
surface t — const of a coordinate system, and consequently is not a priori compatible with the 
Lorentz invariance. Thus we expect that the equations of motion in harmonic coordinates (which 
manifestly preserve the global Lorentz invariance) should exhibit at some stage a violation of the 
Lorentz invariance due to the latter regularization. In fact this occurs exactly at the 3PN order. 
Up to the 2.5PN level, the use of the regularization (F)i is sufficient to get some unambiguous 
equations of motion which are Lorentz invariant 42 . To deal with the problem at 3PN order, a 

23 It was shown in Ref. 1381 that using one or the other of these derivatives results in some equations of motion 
that differ by a mere coordinate transformation. This result indicates that the distributional derivatives introduced 
in Ref. 1361 constitute merely some technical tools which are devoid of physical meaning. 
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Lorentz-invariant variant of the regularization, denoted [F]i, was introduced in Ref. |39j . It consists 
of performing the Hadamard regularization within the spatial hypersurface that is geometrically 
orthogonal (in a Minkowskian sense) to the four- velocity of the particle. The regularization [F]i 
differs from the simpler regularization (F)± by relativistic corrections of order 1/c 2 at least. See 
Ref. for the formulas defining this regularization in the form of some infinite power series in 
1/c 2 . The regularization [F]\ plays a crucial role in obtaining the equations of motion at the 3PN 
order in Refs. |37l 138) . In particular, the use of the Lorentz-invariant regularization [F]\ permits 
to obtain the value of the ambiguity parameter Wkmctic m Equation i]132fl below. 

8.2 Hadamard regularization ambiguities 

The "standard" Hadamard regularization yields some ambiguous results for the computation of 
certain integrals at the 3PN order, as Jaranowski and Schafer |136l 11371 1138] first noticed in 
their computation of the equations of motion within the ADM-Hamiltonian formulation of general 
relativity. By standard Hadamard regularization we mean the regularization based solely on the 
definitions of the partie finie of a singular function, Equation (]122|) . and the partie finie of a 
divergent integral, Equation l]124f) (i.e. without using a theory of pseudo-functions and generalized 
distributional derivatives as proposed in Refs. 36 39 ). It was shown in Refs. |136l 11371 1138) that 
there are two and only two types of ambiguous terms in the 3PN Hamiltonian, which were then 
parametrized by two unknown numerical coefficients astatic and akinetic- 

Motivated by the previous result, Blanchet and Faye |361 131?] introduced their "extended" 
Hadamard regularization, the one we outlined in Section 18.11 This new regularization is mathe- 
matically well-defined and free of ambiguities; in particular it yields unique results for the computa- 
tion of any of the integrals occuring in the 3PN equations of motion. Unfortunately, the extended 
Hadamard regularization turned out to be in a sense incomplete, because it was found |37l 138] 
that the 3PN equations of motion involve one and only one unknown numerical constant, called 
A, which cannot be determined within the method. The comparison of this result with the work 
of Jaranowski and Schafer |136l 1137) . on the basis of the computation of the invariant energy of 
compact binaries moving on circular orbits, showed |37] that 

41 

^kinetic = (132) 

11 x 1987 nrrt 

^static — ~~^~ A ^jjj"- i 1 ^) 

Therefore, the ambiguity Wkinotic is fixed, while A is equivalent to the other ambiguity ecstatic- 
Notice that the value (|132]> for the kinetic ambiguity parameter ^kinetic, which is in factor of 
some velocity dependent terms, is the only one for which the 3PN equations of motion are Lorentz 
invariant. Fixing up this value was possible because the extended Hadamard regularization [3l"l 13*9"] 
was defined in such a way that it keeps the Lorentz invariance. 

Damour, Jaranowski, and Schafer [95] recovered the value of akinetic given in Equation l]132|l 
by directly proving that this value is the unique one for which the global Poincare invariance of 
the ADM-Hamiltonian formalism is verified. Since the coordinate conditions associated with the 
ADM formalism do not manifestly respect the Poincare symmetry, they had to prove that the 3PN 
Hamiltonian is compatible with the existence of generators for the Poincare algebra. By contrast, 
the harmonic-coordinate conditions preserve the Poincare invariance, and therefore the associated 
equations of motion at 3PN order should be manifestly Lorentz-invariant, as was indeed found to 
be the case in Refs. [3*71 155] . 

The appearance of one and only one physical unknown coefficient A in the equations of motion 
constitutes a quite striking fact, that is related specifically with the use of a Hadamard- type reg- 
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ularization 24 . Technically speaking, the presence of the ambiguity parameter A is associated with 
the non-distributivity of Hadamard's regularization, in the sense of Equation i|123fl . Mathemati- 
cally speaking, A is probably related to the fact that it is impossible to construct a distributional 
derivative operator, such as Equations I |129lll30*lll31| l. satisfying the Leibniz rule for the derivation 
of the product j!95j . The Einstein field equations can be written in many different forms, by shift- 
ing the derivatives and operating some terms by parts with the help of the Leibniz rule. All these 
forms are equivalent in the case of regular sources, but since the derivative operator H129lll3()lll31|l 
violates the Leibniz rule they become inequivalent for point particles. Finally, physically speaking, 
let us argue that A has its root in the fact that in a complete computation of the equations of 
motion valid for two regular extended weakly self-gravitating bodies, many non-linear integrals, 
when taken individually ', start depending, from the 3PN order, on the internal structure of the 
bodies, even in the "compact-body" limit where the radii tend to zero. However, when considering 
the full equations of motion, we expect that all the terms depending on the internal structure can 
be removed, in the compact-body limit, by a coordinate transformation (or by some appropriate 
shifts of the central world lines of the bodies), and that finally A is given by a pure number, for 
instance a rational fraction, independent of the details of the internal structure of the compact 
bodies. From this argument (which could be justified by the effacing principle in general relativity) 
the value of A is necessarily the one we compute below, Equation (|135|l . and will be valid for any 
compact objects, for instance black holes. 

The ambiguity parameter w s tatici which is in factor of some static, velocity-independent term, 
and hence cannot be derived by invoking Lorentz invariance, was computed by Damour, Jara- 
nowski, and Schafer |9fi| by means of dimensional regularization, instead of some Hadamard-type 
one, within the ADM-Hamiltonian formalism. Their result is 

^static = 0. (134) 

As Damour et al. 9*5] argue, clearing up the static ambiguity is made possible by the fact that 
dimensional regularization, contrary to Hadamard's regularization, respects all the basic properties 
of the algebraic and differential calculus of ordinary functions: associativity, commutativity and 
distributivity of point-wise addition and multiplication, Leibniz's rule, and the Schwarz lemma. In 
this respect, dimensional regularization is certainly better than Hadamard's one, which does not 
respect the distributivity of the product (recall Equation (|123|^ and unavoidably violates at some 
stage the Leibniz rule for the differentiation of a product. 

The ambiguity parameter A is fixed from the result (|134|) and the necessary link i|133|l provided 
by the equivalence between the harmonic-coordinates and ADM-Hamiltonian formalisms |37l 197) . 
However, A was also been computed directly by Blanchet, Damour, and Esposito-Farese J30J ap- 
plying dimensional regularization to the 3PN equations of motion in harmonic coordinates (in the 
line of Refs. jSZIEHl)- The end result, 

_ 1987 

3080' v ' 

is in full agreement with Equation (|134|) 25 . Besides the independent confirmation of the value of 
^static or A, the work |30| provides also a confirmation of the consistency of dimensional regular- 
ization, because the explicit calculations are entirely different from the ones of Ref. |9fij : Harmonic 



24 Note also that the harmonic-coordinates 3PN equations of motion as they have been obtained in Refs. 1371 1381 
depend, in addition to A, on two arbitrary constants and rL parametrizing some logarithmic terms. These 
constants are closely related to the constants si and S2 in the partie-finie integral 11241 : see Ref. 1381 for the 
precise definition. However, and r' 2 are not "physical" in the sense that they can be removed by a coordinate 
transformation. 

25 One may wonder why the value of A is a complicated rational fraction while ^ stat i c is so simple. This is because 
^static was introduced precisely to measure the amount of ambiguities of certain integrals, while, by contrast, A was 
introduced as an unknown constant entering the relation between the arbitrary scales r[ , r' 2 
si,S2 on the other hand, which has a priori nothing to do with ambiguities of integrals. 
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coordinates are used instead of ADM-type ones, the work is at the level of the equations of motion 
instead of the Hamiltonian, and a different form of Einstein's field equations is solved by a different 
iteration scheme. 

Let us comment here that the use of a self-field regularization, be it dimensional or based 
on Hadamard's partie finie, signals a somewhat unsatisfactory situation on the physical point of 
view, because, ideally, we would like to perform a complete calculation valid for extended bodies, 
taking into account the details of the internal structure of the bodies (energy density, pressure, 
internal velocity field, etc.). By considering the limit where the radii of the objects tend to zero, 
one should recover the same result as obtained by means of the point-mass regularization. This 
would demonstrate the suitability of the regularization. This program was undertaken at the 2PN 
order by Kopcikin et al. [14611124) who derived the equations of motion of two extended fluid balls, 
and obtained equations of motion depending only on the two masses m\ and mi of the compact 
bodies 26 . At the 3PN order we expect that the extended-body program should give the value 
of the regularization parameter A (maybe after some gauge transformation to remove the terms 
depending on the internal structure). Ideally, its value should be confirmed by independent and 
more physical methods (like those of Refs. [2111 11471 TlOl ). 

An important work, in several respects more physical than the formal use of regularizations, is 
the one of Itoh and Futamase [13UII129] . following previous investigations in Refs. [13111132] , These 
authors derived the 3PN equations of motion in harmonic coordinates by means of a particular 
variant of the famous "surface-integral" method introduced long ago by Einstein, Infeld, and 
Hoffmann [106) . The aim is to describe extended relativistic compact binary systems in the strong- 
field point particle limit defined in Ref. [115) . This approach is very interesting because it is based 
on the physical notion of extended compact bodies in general relativity, and is free of the problems 
of ambiguities due to the Hadamard self-field regularization. The end result of Refs. [131)1 1129| is 
in agreement with the 3PN harmonic coordinates equations of motion [371 138) and, moreover, is 
unambiguous, as it does determine the ambiguity parameter A to exactly the value ([135)1 . 

We next consider the problem of the binary's radiation field, where the same phenomenon 
occurs, with the appearance of some Hadamard regularization ambiguity parameters at 3PN or- 
der. More precisely, Blanchet, Iyer, and Joguet , in their computation of the 3PN compact 
binary's mass quadrupole moment Iy, found it necessary to introduce three Hadamard regular- 
ization constants £, k, and £, which are additional to and independent of the equation-of-motion 
related constant A. The total gravitational-wave flux at 3PN order, in the case of circular orbits, 
was found to depend on a single combination of the latter constants, 9 = £ + 2k + C, and the 
binary's orbital phase, for circular orbits, involves only the linear combination of 9 and A given by 
9 = 9 — 7A/3, as shown in |4l)|. 

Dimensional regularization (instead of Hadamard's) has next been applied by Blanchet, Damour, 
Esposito-Farese, and Iyer [311132) to the computation of the 3PN radiation field of compact binaries, 
leading to the following unique values for the ambiguity parameters 27 : 

9871 
^ ~ ~9240' 

« = 0, (136) 



These values represent the end result of dimensional regularization. However, several alterna- 
tive calculations provide a check, independent of dimensional regularization, for all the param- 

26 See some comments on this work in Ref. |84|. pp. 168-169. 

27 The result for £ happens to be amazingly related to the one for A by a cyclic permutation of digits; compare 
3£ = -9871/3080 with A = -1987/3080. 
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eters (|136|) . Blanchct and Iyer |44] compute the 3PN binary's mass dipole moment li using 
Hadamard's regularization, and identify Ij with the 3PN center of mass vector position Gi, al- 
ready known as a conserved integral associated with the Poincare invariance of the 3PN equations 
of motion in harmonic coordinates |103) . This yields £ + n = —9871/9240 in agreement with 
Equation H136JI . Next, we consider [2] the limiting physical situation where the mass of one of 
the particles is exactly zero (say, TO2 = 0), and the other particle moves with uniform velocity. 
Technically, the 3PN quadrupole moment of a boosted Schwarzschild black hole is computed and 
compared with the result for Iy in the limit m,2 = 0. The result is £ = —7/33, and represents a 
direct verification of the global Poincare invariance of the wave generation formalism (the param- 
eter £ represents the analogue for the radiation field of the equation-of-motion related parameter 
^kinetic) 28 - Finally, k = is proven [[""""] by showing that there are no dangerously divergent 
"diagrams" corresponding to non-zero K-values, where a diagram is meant here in the sense of 
Ref. |S7]. 

The determination of the parameters (|136(l completes the problem of the general relativistic 
prediction for the templates of inspiralling compact binaries up to 3PN order (and actually up 
to 3.5PN order as the corresponding tail terms have already been determined The relevant 

combination of the parameters (|136(l entering the 3PN energy flux in the case of circular orbits is 
now fixed to be 

*^ + 2 « + < = -w- < 137 » 

Numerically, 9 ~ —1.28041. The orbital phase of compact binaries, in the adiabatic inspiral regime 
(i.e. evolving by radiation reaction), involves at 3PN order a combination of parameters which is 
determined as 

The fact that the numerical value of this parameter is quite small, 8 ~ 0.22489, indicates, following 
measurement- accuracy analyses [55] "5S] ["^ , that the 3PN (or, even better, 3.5PN) order should 
provide an excellent approximation for both the on-line search and the subsequent off-line analysis 
of gravitational wave signals from inspiralling compact binaries in the LIGO and VIRGO detectors. 



8.3 Dimensional regularization of the equations of motion 

As reviewed in Section work at 3PN order using Hadamard's self-field regularization showed 
the appearance of ambiguity parameters, due to an incompleteness of the Hadamard regularization 
employed for curing the infinite self field of point particles. We give here more details on the 
determination using dimensional regularization of the ambiguity parameter A which appeared in 
the 3PN equations of motion (recall that A is equivalent to the static ambiguity parameter w s t a tic, 
see Equation (|133(l ). 

Dimensional regularization was invented as a means to preserve the gauge symmetry of per- 
turbative quantum field theories |199l 1511 1571 173*] . Our basic problem here is to respect the gauge 
symmetry associated with the diffcomorphism invariance of the classical general relativistic de- 
scription of interacting point masses. Hence, we use dimensional regularization not merely as a 
trick to compute some particular integrals which would otherwise be divergent, but as a powerful 
tool for solving in a consistent way the Einstein field equations with singular point-mass sources, 
while preserving its crucial symmetries. In particular, we shall prove that dimensional regulariza- 
tion determines the kinetic ambiguity parameter akinetic (and its radiation-field analogue C) , and is 

28 The work I34| provided also some new expressions for the multipole moments of an isolated post-Newtonian 
source, alternative to those given by Theorem"""] in the form of surface integrals extending on the outer part of the 
source's near zone. 
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therefore able to correctly keep track of the global Lorentz-Poincare invariance of the gravitational 
field of isolated systems. 

The Einstein field equations in d+1 space-time dimensions, relaxed by the condition of harmonic 
coordinates d^h a ^ = 0, take exactly the same form as given in Equations 114|) . In particular □ 
denotes the flat space-time d'Alembertian operator in d-l- 1 dimensions. The gravitational constant 
G is related to the usual three-dimensional Newton's constant Gat by 

G = G N £%- 3 , (139) 

where £o denotes an arbitrary length scale. The explicit expression of the gravitational source term 
A a P involves some d-dependent coefficients, and is given by 

- g^gurdxh^d^ - g f, »g VT d x h ar d lt h v * + g^g XT d x h^d T h^ 

+ \{2g a »g^ - g^gn [g XT g^ - jL^g^g^ d^d v h T \ (140) 

When d = 3 we recover Equation l|15f) . In the following we assume, as usual in dimensional 
regularization, that the dimension of space is a complex number, d € C, and prove many results 
by invoking complex analytic continuation in d. We shall pose s = d — 3. 

We parametrize the 3PN metric in d dimensions by means of straightforward d-dimensional 
generalizations of the retarded potentials V, Vi, Wij, Ri, and X of Section [7| Those are obtained 
by post-Newtonian iteration of the d-dimensional field equations, starting from the following defi- 
nitions of matter source densities 

2 (d - 2)T 00 + T u 

a — 



d-l c 2 

T 0i (141) 

0i = , 

c 

which generalize Equations l|116Jl . As a result all the expressions of Section acquire some explicit 
(i-dependent coefficients. For instance we find [23 

V = n-\ l-inGa] , 



\ i (142) 



Here means the retarded integral in d + 1 space-time dimensions, which admits, though, no 
simple expression in physical (t, x) space. 

As reviewed in Section f8. II the generic functions we have to deal with in 3 dimensions, say F(x), 
are smooth on M 3 except at yi and y2 , around which they admit singular Laurent-type expansions 
in powers and inverse powers of r\ = |x — yi| and = |x — y2 1 , given by Equation l|121ll . In d 
spatial dimensions, there is an analogue of the function F, which results from the post-Newtonian 
iteration process performed in d dimensions as we just outlined. Let us call this function F^ d \x) 7 
where x G WL d . When r\ — > the function F^ d ' admits a singular expansion which is a little bit 
more complicated than in 3 dimensions, as it reads 

FW(x) = J2 r^fi^i) + o(»f )■ (143) 

Po<P<N 1 

(/()<<-/<<-/! 
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The coefficients f x p^g(ni) depend on e = d — 3, and the powers of r\ involve the relative integers 
p and q whose values are limited by some po, qo, and q\ as indicated. Here we will be interested 
in functions F^ d \x) which have no poles as e — » (this will always be the case at 3PN order). 
Therefore, we can deduce from the fact that p( d )(x) is continuous at d = 3 the constraint 

£/^ =0) (ni) = / P (ni). (144) 

9=90 

For the problem at hand, we essentially have to deal with the regularization of Poisson integrals, 
or iterated Poisson integrals (and their gradients needed in the equations of motion) , of the generic 
function fw. The Poisson integral of F^ d \ in d dimensions, is given by the Green's function for 
the Laplace operator, 

pW(x') = A- 1 |>«> (X )1 =-A / dd * fW (x) , (145) 



4ir _ 

where k is a constant related to the usual Eulerian T-function by 29 



k=-W-. (146) 

TT 2 

We need to evaluate the Poisson integral at the point x' = yi where it is singular; this is quite easy 
in dimensional regularization, because the nice properties of analytic continuation allow simply to 
get [p( d ) (x')] x / =yi by replacing x' by yi in the explicit integral form (|145|) . So we simply have 

^ rf) (yx) = -^/^ (d) (x). (147) 

It is not possible at present to compute the equations of motion in the general <i-dimensional 
case, but only in the limit where e — > [HHHS!]- The main technical step of our strategy consists of 
computing, in the limit £ — > 0, the difference between the d-dimensional Poisson potential (|147|) . 
and its Hadamard 3-dimensional counterpart given by (P)i, where the Hadamard partie finie 
is defined by Equation (|122|) . Actually, we must be very precise when defining the Hadamard 
partie finie of a Poisson integral. Indeed, the definition 1|122(1 stricto sensu is applicable when the 
expansion of the function F, when r\ — > 0, does not involve logarithms of r\\ see Equation (|121|) . 
However, the Poisson integral P(x') of F(x) will typically involve such logarithms at the 3PN order, 
namely some lnr^ where r' x = |x' — yi| formally tends to zero (hence lnr^ is formally infinite). The 
proper way to define the Hadamard partie finie in this case is to include the In r' x into its definition, 
so we arrive at [3^] 

( p )i = ~ Pfri, S2 J ^F(x) {r\ F) x . (148) 

The first term follows from Hadamard's partie finie integral (|124|) : the second one is given by 
Equation (|122(l . Notice that in this result the constant si entering the partie finie integral Ijl24|l 
has been "replaced" by r[, which plays the role of a new regularization constant (together with 
r' 2 for the other particle), and which ultimately parametrizes the final Hadamard regularized 3PN 

29 We have lim t ;_,3 k = 1. Notice that k is closely linked to the volume fld—l °f the sphere with d — 1 dimensions 
(i.e. embedded into Euclidean d-dimensional space): 

47T 

kil d _, = -. 
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equations of motion. It was shown that r[ and r' 2 are unphysical, in the sense that they can be 
removed by a coordinate transformation On the other hand, the constant s 2 remaining in 

the result i|148fl is the source for the appearance of the physical ambiguity parameter A, as it will be 
related to it by Equation (|150() . Denoting the difference between the dimensional and Hadamard 
regularizations by means of the script letter T>, we pose (for the result concerning the point 1) 

PP 1 =pW(yi)-(P) 1 . (149) 

That is, T>P\ is what we shall have to add to the Hadamard-regularization result in order to get the 
d-dimensional result. However, we shall only compute the first two terms of the Laurent expansion 
of T)P\ when e — ► 0, say a_ie _1 + ao + 0(e). This is the information we need to clear up the 
ambiguity parameter. We insist that the difference T> P\ comes exclusively from the contribution 
of terms developing some poles oc 1/e in the <i-dimensional calculation. 

Next we outline the way we obtain, starting from the computation of the "difference" , the 3PN 
equations of motion in dimensional regularization, and show how the ambiguity parameter A is 
determined. By contrast to r[ and r 2 which are pure gauge, A is a genuine physical ambiguity, 
introduced in Refs. |36l I38| as the single unknown numerical constant parametrizing the ratio 
between s 2 and r' 2 (where s 2 is the constant left in Equation 1148|l ) as 

f r' 2 \ 159 m 1 +m 2 

ln — =^ + A and 1^2, 150 

\s 2 J 308 m 2 

where mi and m 2 are the two masses. The terms corresponding to the A-ambiguity in the accel- 
eration ai = dvi/dt of particle 1 read simply 

44A G% mi ml (mi + m 2 ) , . 

Aai [A] = -g— g m 2 , (151) 

O 1 12 C 

where the relative distance between particles is denoted yi — y 2 = r\ 2 ni2 (with ni2 being the unit 
vector pointing from particle 2 to particle 1). We start from the end result of Ref. |3H1 for the 3PN 
harmonic coordinates acceleration ai in Hadamard's regularization, abbreviated as HR. Since the 
result was obtained by means of the specific extended variant of Hadamard's regularization (in 
short EHR, see Section ISTTft we write it as 

a (HR) =a (EHR) +Aai[A]) (152) 

where ap HR '' is a fully determined functional of the masses mi and m-2, the relative distance 
T\2^i2, the coordinate velocities vi and V2, and also the gauge constants r[ and r 2 . The only 
ambiguous term is the second one and is given by Equation 115111 . 

Our strategy is to express both the dimensional and Hadamard regularizations in terms of 
their common "core" part, obtained by applying the so-called "pure-Hadamard-Schwartz" (pHS) 
regularization. Following the definition of Ref. |3Uj . the pHS regularization is a specific, minimal 
Hadamard- type regularization of integrals, based on the partie finie integral i|124fl . together with a 
minimal treatment of "contact" terms, in which the definition i|124fl is applied separately to each 
of the elementary potentials V, Vi, . . . (and gradients) that enter the post-Newtonian metric in the 
form given in Section[7] Furthermore, the regularization of a product of these potentials is assumed 
to be distributive, i.e. (FG)i — (F)i(G)i in the case where F and G are given by such elementary 
potentials (this is in contrast with Equation (|123|l ). The pHS regularization also assumes the use 
of standard Schwartz distributional derivatives |196) . The interest of the pHS regularization is that 
the dimensional regularization is equal to it plus the "difference"; see Equation l|155fl . 

To obtain the pHS-regularized acceleration we need to substract from the EHR result a series 
of contributions, which are specific consequences of the use of EHR [3(312111. For instance, one 
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of these contributions corresponds to the fact that in the EHR the distributional derivative is 
given by Equations i)129l I130|) which differs from the Schwartz distributional derivative in the pHS 
regularization. Hence we define 

^PHS^EHR)^^ (153) 
A 

where the Sa&iS denote the extra terms following from the EHR prescriptions. The pHS-regularized 
acceleration (1153)1 constitutes essentially the result of the first stage of the calculation of ai, as 
reported in Ref. |1U9| . 

The next step consists of evaluating the Laurent expansion, in powers of e = d — 3, of the 
difference between the dimensional regularization and the pHS (3-dimensional) computation. As 
we reviewed above, this difference makes a contribution only when a term generates a pole ~ 1/e, 
in which case the dimensional regularization adds an extra contribution, made of the pole and the 
finite part associated with the pole (we consistently neglect all terms O(e)). One must then be 
especially wary of combinations of terms whose pole parts finally cancel ( "cancelled poles" ) but 
whose dimensionally regularized finite parts generally do not, and must be evaluated with care. 
We denote the above defined difference by 

Va 1 =^2vP 1 . (154) 

It is made of the sum of all the individual differences of Poisson or Poisson-like integrals as com- 
puted in Equation (|149|l . The total difference (|154|) depends on the Hadamard regularization scales 
r[ and S2 (or equivalently on A and r[, r' 2 ), and on the parameters associated with dimensional reg- 
ularization, namely e and the characteristic length scale £q introduced in Equation (|139() . Finally, 
our main result is the explicit computation of the e-expansion of the dimensional regularization 
(DR) acceleration as 

a (DR) =a (pHS) +pa ^ (155) 

With this result we can prove two theorems |3Uj: 

Theorem 8 The pole part oc 1/e of the DR acceleration I J. 5.5)) can be re-absorbed (i.e. renormal- 
ized) into some shifts of the two "bare" world-lines: yi — * yi + £i and y2 — * y2 + £2? with, say, 
€1.2 oc 1/e, so that the result, expressed in terms of the "dressed" quantities, is finite when e — ► 0. 

The situation in harmonic coordinates is to be contrasted with the calculation in ADM-type coor- 
dinates within the Hamiltonian formalism, where it was shown that all pole parts directly cancel 
out in the total 3PN Hamiltonian: No renormalization of the world-lines is needed • A central 
result is then as follows: 



Theorem 9 The renormalized (finite) DR acceleration is physically equivalent to the Hadamard- 
regularized (HR) acceleration (end result of Ref. (381), in the sense that 



a (HR) _ 



(DR) , c 

a\ ' + d£ ai 



(156) 



where 6^ sli denotes the effect of the shifts on the acceleration, if and only if the HR ambiguity 
parameter X entering the harmonic- coordinates equations of motion takes the unique value 1) 135\) . 



The precise shifts £1 and £2 needed in TheoremElinvolve not only a pole contribution oc 1/e (which 
would define a renormalization by minimal subtraction (MS)), but also a finite contribution when 
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0. Their explicit expressions read 30 : 
11 G? 



€1 = 




327 
1540 



ajvi (together with 1 «-> 2), 



(157) 



where Gn is Newton's constant, £o is the characteristic length scale of dimensional regularization 
(cf. Equation l|139[) ). ajvi is the Newtonian acceleration of the particle 1 in d dimensions, and 
q = 4ire c depends on Euler's constant C = 0.577 



8.4 Dimensional regularization of the radiation field 

We now address the similar problem concerning the binary's radiation field (3PN beyond the 
Einstein quadrupole formalism), for which three ambiguity parameters, £, k, £, have been shown 
to appear 45, HI] (see Section 

To apply dimensional regularization, we must use as in Section 18.31 the d-dimensional post- 
Newtonian iteration [leading to equations such as l|142|l ] ; and, crucially, we have to generalize to d 
dimensions some key results of the wave generation formalism of Part [I] Essentially we need the 
d-dimensional analogues of the multipole moments of an isolated source II and Jl, Equations H85JI . 
The result we find in the case of the mass-type moments is 

2(d + 2l-2) A (2 ) 1 

+ c^d + l-l)(d + l-2)(d + 2l + 2) XabL { ^ 2 f (X ' t] J ' (158) 

where we denote (generalizing Equations JHEJ) 

2 (d - 2)r u 



£ = 



1 



(159) 



c 



13 



and where for any source densities the underscript [I] means the infinite series 

E(x , t ) = y * Lii±*) , ( M o) 2 s( X) t) . (i6o) 

The latter definition represents the d-dimensional version of the post-Newtonian expansion se- 
ries H9H . At Newtonian order, Equation (|158(l reduces to the standard result = J d d x.pxL + 
0{cr 2 ) withp = T 00 /c 2 . 

The ambiguity parameters £, k, and £ come from the Hadamard regularization of the mass 
quadrupole moment Iy at the 3PN order. The terms corresponding to these ambiguities were 
found to be 



«,C] 



1 <-> 2, (161) 



mi 

3(> Whcn working at the level of the equations of motion (not considering the metric outside the world-lines), the 
effect of shifts can be seen as being induced by a coordinate transformation of the bulk metric as in Ref. |38| . 
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where yi, vi, and ai denote the first particle's position, velocity, and acceleration. We recall that 
the brackets ( ) surrounding indices refer to the symmetric-trace-free (STF) projection. Like in 
Section [8.3l we express both the Hadamard and dimensional results in terms of the more basic pHS 
regularization. The first step of the calculation 01] is therefore to relate the Hadamard-regularized 
quadrupole moment llf R \ for general orbits, to its pHS part: 



i j 

(HR) = (pHS) M 



\ 1 , 9 



(162) 



In the right-hand side we find both the pHS part, and the effect of adding the ambiguities, with 
some numerical shifts of the ambiguity parameters coming from the difference between the specific 
Hadamard-type regularization scheme used in Ref. and the pHS one. The pHS part is free of 
ambiguities but depends on the gauge constants r[ and r' 2 introduced in the harmonic-coordinates 
equations of motion |37l 138) . 

We next use the d-dimensional moment (|158|) to compute the difference between the dimensional 
regularization (DR) result and the pHS one |31l I32| . As in the work on equations of motion, we 
find that the ambiguities arise solely from the terms in the integration regions near the particles 
(i.e. r\ = |x — yi| — > or r 2 = |x — y2 1 — > 0) that give rise to poles cx 1/e, corresponding to 
logarithmic ultra-violet (UV) divergences in 3 dimensions. The infra-red (IR) region at infinity 
(i.e. |x| — > +oo) does not contribute to the difference DR — pHS. The compact-support terms in 
the integrand of Equation (|158|) . proportional to the matter source densities a, o~ a , and o~ a b, are 
also found not to contribute to the difference. We are therefore left with evaluating the difference 
linked with the computation of the non-compact terms in the expansion of the integrand in l|158f) 
near the singularities that produce poles in d dimensions. 

Let F^(x) be the non-compact part of the integrand of the quadrupole moment l|158[) (with 
indices L = ij), where F^ includes the appropriate multipolar factors such cis Xij j so that 



1^ 



J Af'^x). (163) 



We do not indicate that we are considering here only the non-compact part of the moments. Near 
the singularities the function F^ d \x) admits a singular expansion of the type H143[) . In practice, 
the various coefficients i/pfg are computed by specializing the general expressions of the non- linear 
retarded potentials V, V a , W a b, ■ ■ ■ (valid for general extended sources) to the point particles case 
in d dimensions. On the other hand, the analogue of Equation l|163f) in 3 dimensions is 



= Pf J d 3 xF(x), (164) 



where Pf refers to the Hadamard partie finic defined in Equation (|124|) . The difference T>I between 
the DR evaluation of the d-dimensional integral (|163|) . and its corresponding three-dimensional 
evaluation, i.e. the partie finic 1164fl . reads then 

i,,. (165) 

Such difference depends only on the UV behaviour of the integrands, and can therefore be computed 
"locally" , i.e. in the vicinity of the particles, when r\ — > and r 2 — > 0. We find that Equation l|165fl 
depends on two constant scales s\ and S2 coming from Hadamard's partie finie l|124|l . and on the 
constants belonging to dimensional regularization, which are e — d — 3 and the length scale £q 
defined by Equation (|139(l . The dimensional regularization of the 3PN quadrupole moment is then 
obtained as the sum of the pHS part, and of the difference computed according to Equation (|165|l . 
namely 

I (DR )=I (pHS) +pI ^ (16g) 
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An important fact, hidden in our too-compact notation ljl66|) . is that the sum of the two terms in 
the right-hand side of Equation i|166[l does not depend on the Hadamard regularization scales Si 
and S2- Therefore it is possible without changing the sum to re-express these two terms (separately) 
by means of the constants r[ and r' 2 instead of si and S2, where r[, r' 2 are the two fiducial scales 
entering the Hadamard-regularization result H162J1 • This replacement being made the pHS term in 
Equation l|166[l is exactly the same as the one in Equation (|162|) . At this stage all elements are in 
place to prove the following theorem 02] : 

Theorem 10 The DR quadrupole moment is physically equivalent to the Hadamard- 

regularized one (end result of Refs. \4o\ \44\), in the sense that 



ig^ = Hm 



(167) 



where 5^1^ denotes the effect of the same shifts as determined in Theorems^ and\Q if and only 
if the HR ambiguity parameters k, and £ take the unique values \13b]) . Moreover, the poles 
1/e separately present in the two terms in the brackets of Equation J_?6'7| ) cancel out, so that the 
physical ( "dressed" ) DR quadrupole moment is finite and given by the limit when e — > as shown 
in Equation \ lbty . 

This theorem finally provides an unambiguous determination of the 3PN radiation field by dimen- 
sional regularization. Furthermore, as reviewed in Section 18.21 several checks of this calculation 
could be done, which provide, together with comparisons with alternative methods |96l 1301 11301 
I129J . independent confirmations for the four ambiguity parameters A, £, k, and and confirm 
the consistency of dimensional regularization and its validity for describing the general-relativistic 
dynamics of compact bodies. 
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9 Newtonian-like Equations of Motion 



9.1 The 3PN acceleration and energy 

We present the acceleration of one of the particles, say the particle 1, at the 3PN order, as well 
as the 3PN energy of the binary, which is conserved in the absence of radiation reaction. To get 
this result we used essentially a "direct" post-Newtonian method (issued from Ref. 03), which 
consists of reducing the 3PN metric of an extended regular source, worked out in Equations (jll5|l . 
to the case where the matter tensor is made of delta functions, and then curing the self-field 
divergences by means of the Hadamard regularization technique. The equations of motion are 
simply the geodesic equations associated with the regularized metric (see Ref. for a proof). 
The Hadamard ambiguity parameter A is computed from dimensional regularization in Section f8.3l 
We also add the 3.5PN terms which are known from Refs. (1531 IHH IHBl IT7T1 IH51 ITHT] . 

Though the successive post-Newtonian approximations are really a consequence of general 
relativity, the final equations of motion must be interpreted in a Newtonian-like fashion. That is, 
once a convenient general-relativistic (Cartesian) coordinate system is chosen, we should express 
the results in terms of the coordinate positions, velocities, and accelerations of the bodies, and 
view the trajectories of the particles as taking place in the absolute Euclidean space of Newton. 
But because the equations of motion are actually relativistic, they must 

(i) stay manifestly invariant - at least in harmonic coordinates - when we perform a global 
post-Newtonian-expanded Lorentz transformation, 

(ii) possess the correct "perturbative" limit, given by the geodesies of the (post-Newtonian- 
expanded) Schwarzschild metric, when one of the masses tends to zero, and 

(hi) be conservative, i.e. to admit a Lagrangian or Hamiltonian formulation, when the gravita- 
tional radiation reaction is turned off. 

We denote by ri 2 = |yi(£) — y2 (*) | the harmonic-coordinate distance between the two particles, 
with yi = (y\) and y2 = {1)2)1 by n\ 2 — (y[ — y 2 )/ri 2 the corresponding unit direction, and by 
v\ = dy\/dt and a\ = dv\/dt the coordinate velocity and acceleration of the particle 1 (and idem 
for 2). Sometimes we pose v\ 2 — v\ — v\ for the relative velocity. The usual Euclidean scalar 
product of vectors is denoted with parentheses, e.g., (ni 2 «i) = ni 2 • vi and (vi« 2 ) = Vi • v 2 . The 
equations of the body 2 are obtained by exchanging all the particle labels 1^2 (remembering 
that n\ 2 and v\ 2 change sign in this operation): 
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The 2.5PN and 3.5PN terms are associated with gravitational radiation reaction. The 3PN 
harmonic-coordinates equations of motion depend on two arbitrary length scales r[ and r 2 as- 
sociated with the logarithms present at the 3PN order 31 . It has been proved in Ref. that r[ 
and r' 2 are merely linked with the choice of coordinates - we can refer to r[ and r 2 as "gauge 
constants". In our approach |37l 138). the harmonic coordinate system is not uniquely fixed by 
the coordinate condition d^h 01 ^ = 0. In fact there are infinitely many harmonic coordinate sys- 
tems that are local. For general smooth sources, as in the general formalism of Part[l] we expect 

31 Notice also the dependence upon it 2 . Technically, the it 2 terms arise from non-linear interactions involving some 
integrals such as 



1 



ri2 
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the existence and uniqueness of a global harmonic coordinate system. But here we have some 
point-particles, with delta-function singularities, and in this case we do not have the notion of a 
global coordinate system. We can always change the harmonic coordinates by means of the gauge 
vector r] a — Sx a , satisfying A?/ Q = except at the location of the two particles (we assume that 
the transformation is at the 3PN level, so we can consider simply a flat-space Laplace equation). 
More precisely, we can show that the logarithms appearing in Equation (|168|1 . together with the 
constants r[ and r' 2 therein, can be removed by the coordinate transformation associated with the 
3PN gauge vector (with n = |x — yi(t)| and r 2 = |x — y 2 (i)|): 
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(169) 

Therefore, the "ambiguity" in the choice of the constants r[ and r' 2 is completely innocuous on 
the physical point of view, because the physical results must be gauge invariant. Indeed we shall 
verify that r[ and r' 2 cancel out in our final results. 

When retaining the "even" relativistic corrections at the 1PN, 2PN and 3PN orders, and 
neglecting the "odd" radiation reaction terms at the 2.5PN and 3.5PN orders, we find that the 
equations of motion admit a conserved energy (and a Lagrangian, as we shall see), and that energy 
can be straightforwardly obtained by guess-work starting from Equation (|168|) , with the result 
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To the terms given above, we must add the terms corresponding to the relabelling 1^2. Actually, 
this energy is not conserved because of the radiation reaction. Thus its time derivative, as computed 
by means of the 3PN equations of motion themselves (i.e. order-reducing all the accelerations), is 
purely equal to the 2.5PN effect, 



dE _ 4 G z mim 2 
dt 5 c 5 rf 2 



(vivi-i 



t (-in 1 1 ^ ( Vm-2^ 



r\ 2 



rv2 



{ni 2 vi)(n 12 vi 2 ) [3v 



,,2 



Gm-t 52 Gm,2 
6 - + 



?"12 



3 ri2 
(171 ! 



The resulting "balance equation" can be better expressed by transfering to the left-hand side 
certain 2.5PN terms so that the right-hand side takes the familiar form of a total energy flux. 
Posing 
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we find agreement with the standard Einstein quadrupole formula Q |SJ) : 
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where the Newtonian trace-free quadrupole moment is Qy = mi(y\y J 1 — hS^yf) + 1 •*->• 2. We refer 
to Iyer and Will |133l I134| for the discussion of the energy balance equation at the next 3.5PN 
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order. As we can see, the 3.5PN equations of motion l|168|) are highly relativistic when describing 
the motion, but concerning the radiation they are in fact 1PN, because they contain merely the 
radiation reaction force at the 2.5PN + 3.5PN orders. 



9.2 Lagrangian and Hamiltonian formulations 

The conservative part of the equations of motion in harmonic coordinates (|168|) is derivable from a 
generalized Lagrangian, depending not only on the positions and velocities of the bodies, but also 
on their accelerations: a\ = dv\/dt and a\ = dv\/dt. As shown by Damour and Deruelle [EE]) the 
accelerations in the harmonic-coordinates Lagrangian occur already from the 2PN order. This fact 
is in accordance with a general result of Martin and Sanz |155j that iV-body equations of motion 
cannot be derived from an ordinary Lagrangian beyond the 1PN level, provided that the gauge 
conditions preserve the Lorentz invariance. Note that we can always arrange for the dependence of 
the Lagrangian upon the accelerations to be linear, at the price of adding some so-called "multi- 
zero" terms to the Lagrangian, which do not modify the equations of motion (see, e.g., Ref. |98|L 
At the 3PN level, we find that the Lagrangian also depends on accelerations. It is notable that 
these accelerations are sufficient - there is no need to include derivatives of accelerations. Note also 
that the Lagrangian is not unique because we can always add to it a total time derivative dF/dt, 
where F depends on the positions and velocities, without changing the dynamics. We find jl()3| 



L 



2r 12 1 2 






C G 2 m\m 2 


m-iVi Gmim-2 
+ -^ + l —l 


c 2 ' 




8 r 12 




G 3 m\m 2 


l9G 3 mlml 


M 


2r? 2 + 


8r? 2 



1 3 7 



G 2 m\m 2 (1 2 7. .. . 1. 2 1 2 7, s 7 2 

-(m 2 wi) - -{n\ 2 vx){nx 2 v 2 ) + -(n 12 v 2 ) + -v x - 7 (uiu 2 ) + t" 2 



r\ 2 \2 X '"^~ L > 2 v ~ x *~ lJV ~ x *~" ' 2 K " l *~* J ' 4" 1 4^-"-' 1 4' 

Gm\m 2 ( 3 . , 2/ \2 7. , 2 i 7 4 3 . . . . . . 

77t(«12«i) {n\ 2 v 2 ) ~ -{n 12 v 2 ) v x + -v x + -{ni 2 vi)(ni 2 v 2 )(viv 2 ) 

7*12 \ 16 8 8 4 

-2v 2 (viv 2 ) + l(v lV2 ) 2 + ^vlvi] + "' J ' 1 



16 1 V 16 

7 1 7 

+ Gm 1 m 2 [--{aiv 2 ){ni 2 v 2 ) - -(ni 2 ai)(ni 2 v 2 ) 2 + -{ni 2 ai)v 2 

1 J G 2 m\m 2 ( 13 4 83 3 35 2 2 245 2 2 

lo lo ^4 

179. . \2 235. \ 22 373 4 529 \ 2 / \ 

+ -j2-(ni2Ui)(rn2«2)Ui - -^-(n 12 v 2 ) v x + —v r + — (m 2 «ij (ui«2) 



97, w \ / , 719 2 , , 463, , 2 7 . 

- y(ni 2 ui)(ni 2 « 2 )(uiu 2 ) - — v^viv^ + —{viv 2 ) - ^(niaui 

1 , X/ s 2 1/ ^2 2 463 2 2 19/ x 2 45 

- 2( ni2W i)( n i2 u 2)w 2 + j(ni 2 v 2 ) u 2 + -^r v i v 2 ~ -^{VlV2)V 2 + — ' 

5mi 

H — 

128 

/3 5 1 

+ GTOim 2 ( -(a 1 u 2 )(n 12 u 1 )(m 2 u 2 ) 2 + — {axv^n^v^ 3 + -{n^a^n^v^in^v^,) 3 



63 



+ 7^("-i2ai)(«i2W2) 4 + ^r(a 1 v 1 )(n 12 v 2 )vl - {a 1 v 2 )(n 1 2V 2 )vf 
16 4 

-2(a 1 v 1 )(ni2V2)(viv 2 ) + ^(aiu 2 )(ni 2 t>2)(«iU2) 

3 5 15 

+ g(«i2ai)(n 12 t>2) 2 (wiW2) - g(^i2ai)("-i2Wi) 2 w 2 + y (ai^i)(«i2^2)^2 

15 1 

— ^(aiv 2 )(n 12 v 2 )v 2 - -(n 12 a 1 )(n 12 w 1 )(n 12 w 2 )w| 



5 

1G 



(n 12 ai)(ni2«2) 2 V2 



G 2 m 2 m 2 { 235, w , 29 , w ,, 235, 
+ — — I - -^j-(a 2 vi)(ni 2 wi) - 24(ni 2 a 2 )(ni2Ui) - -^-(aiv 2 )(n 12 v 2 ) 

17. . , ,o 185 , . n 235 , , 9 

r-(ni2ai)(ni2«2j + -r?r{ni2ai)v 1 —(ni 2 a 2 )v 1 

6 Id 48 

185, . 20, , 2 

—{ni 2 ai)(viv 2 ) + —{n 12 ai)v 2 

8 o 

+ Gmim2 f- A( ni2Wl ) 3 ( ni2V2 ) 3 + i(ni 2 wi)(ni 2 w 2 ) 3 w 2 + ^(n 12 v 2 ) 4 vl 
r\ 2 \ 61 8 8 

-^{n 1 2V 1 )(n 1 2V 2 )v i + 7(«12«2) 2 Wi + -Ui 
16 4 16 

15 

-■^(ni2Wi) 2 («i2W 2 ) 2 (wiW 2 ) + (n-i2fi)(n-i2W2)w 2 (wiW2) 
3 13 5 

+ o( n 12l , 2) 2 ^ 2 0lf 2 ) - TXVx{v\V2) + 777(«l 2 Wl)(ni 2 W 2 )(uiW 2 ) 2 

8 16 16 

+ ^{v!V 2 ) 3 - ^(n 12 V!) 2 v 2 vj - —{n 12Vl ){n 12 v 2 )v 2 v 2 + —v\ 
1 



-g^ 2 («l« 2 )uf 

3G 4 mjm 2 G 4 m\m 2 2 ( 9707 22, /r 12 

1 — In — t 



8rf 2 r\ 2 \ 420 

G 3 m?ml /383, x2 889, , 123, , 2 305 , 

+ 13 ( ^-("12^1) - ^g-(^12Ul)(ni2W 2 ) - — (ni2Wi)(ni2«l 2 )7T 



12 



41 2 , ,439, 

^TT ( UlUl2 ) + _ ( Vl V 2 , 




,2 



G 3 m\m 2 ( 8243, , 2 15541 , ,, .3, n2 15611 
+ Is 1~ -^T7ri n ^vi) + ^^-{ni2Vi){n 12 v 2 ) + -{n 12 v 2 ) + -J^q'^i 

22(n 12 vi)(n 12 v 12 ) In ( 



+1^2 + ofij. (174) 

Witness the accelerations occuring at the 2PN and 3PN orders; see also the gauge-dependent 
logarithms of ri 2 /r[ and r\ 2 jr' 2 . We refer to |1()3) for the explicit expressions of the ten conserved 
quantities corresponding to the integrals of energy (also given in Equation <|170J) ). linear and 
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angular momenta, and center-of-mass position. Notice that while it is strictly forbidden to replace 
the accelerations by the equations of motion in the Lagrangian, this can and should be done in the 
final expressions of the conserved integrals derived from that Lagrangian. 

Now we want to exhibit a transformation of the particles dynamical variables - or contact 
transformation, as it is called in the jargon - which transforms the 3PN harmonic-coordinates 
Lagrangian (|1 T4|> into a new Lagrangian, valid in some ADM or ADM-like coordinate system, and 
such that the associated Hamiltonian coincides with the 3PN Hamiltonian that has been obtained 
by Damour, Jaranowski, and Schafer In ADM coordinates the Lagrangian will be "ordinary", 
depending only on the positions and velocities of the bodies. Let this contact transformation be 
Yi(i) = y\{t) + Sy\(t) and 1^2, where Y± and y\ denote the trajectories in ADM and harmonic 
coordinates, respectively. For this transformation to be able to remove all the accelerations in the 
initial Lagrangian £ harm up to the 3PN order, we determine |103| it to be necessarily of the form 



Sy{ = — 

mi 



<9 L harm OF 1 i 



(and idem 1 <-> 2), where F is a freely adjustable function of the positions and velocities, made of 
2PN and 3PN terms, and where X\ represents a special correction term, that is purely of order 
3PN. The point is that once the function F is specified there is a unique determination of the 
correction term X{ for the contact transformation to work (see Ref. |1(J3| for the details). Thus, 
the freedom we have is entirely coded into the function F, and the work then consists in showing 
that there exists a unique choice of F for which our Lagrangian L harm is physically equivalent, via 
the contact transformation (|175|1 . to the ADM Hamiltonian of Ref. [HH]. An interesting point is 
that not only the transformation must remove all the accelerations in L harm , but it should also 
cancel out all the logarithms ln(ri2/r' 1 ) and ln(ri 2 /r 2 ), because there are no logarithms in ADM 
coordinates. The result we find, which can be checked to be in full agreement with the expression 
of the gauge vector in Equation l|169(l , is that F involves the logarithmic terms 
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(176) 



together with many other non-logarithmic terms (indicated by dots) that are entirely specified by 
the isometry of the harmonic and ADM descriptions of the motion. For this particular choice of 
F the ADM Lagrangian reads 

X r harm r r harm jp / 1 \ 

Inserting into this equation all our explicit expressions we find 
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The notation is the same as in Equation (1 174(1 , except that we use upper-case letters to denote the 
ADM-coordinates positions and velocities; thus, for instance Ni 2 = (Yi — Y 2 )/Ri 2 and (Ni 2 Vi) — 
Ni 2 • Vi. The Hamiltonian is simply deduced from the latter Lagrangian by applying the usual 
Legendre transformation. Posing P{ = dL ADM /dVJ and 1 <-> 2, we get PbllHTllT^l 193111173 ] 3 2 
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Arguably, the results given by the ADM-Hamiltonian formalism (for the problem at hand) look 
simpler than their harmonic-coordinate counterparts. Indeed, the ADM Lagrangian is ordinary - 
no accelerations - and there are no logarithms nor associated gauge constants r[ and r' 2 . But of 
course, one is free to describe the binary motion in whatever coordinates one likes, and the two 
formalisms, harmonic Ijl74|l and ADM (|178l I179")l . describe rigorously the same physics. On the 
other hand, the higher complexity of the harmonic-coordinates Lagrangian (1 174(1 enables one to 
perform more tests of the computations, notably by inquiring about the future of the constants r[ 
and r 2 , that we know must disappear from physical quantities such as the center-of-mass energy 
and the total gravitational-wave flux. 



9.3 Equations of motion in the center-of-mass frame 

In this section we translate the origin of coordinates to the binary's center-of-mass by imposing 
that the binary's dipole I; = (notation of Part^J. Actually the dipole moment is computed as 
the center-of-mass conserved integral associated with the boost symmetry of the 3PN equations of 
motion and Lagrangian |1(J3I I43| . This condition results in the (3PN-accurate, say) relationship 
between the individual positions in the center-of-mass frame y\ and y 2 , and the relative position 
x l = y\ — y\ and velocity v % = v\ — v 2 — dx 1 jdt (formerly denoted y\ 2 and v\ 2 ). We shall also use 
the orbital separation r = |x|, together with n = x/r and f = n ■ v. Mass parameters are the total 
mass to = iiij + rn 2 (to = M in the notation of Part [IJ , the mass difference 8m = mi — m 2 , the 
reduced mass // = to 1 to 2 /to, and the very useful symmetric mass ratio 



/i _ m\m 2 



to (toi + m 2 ) 2 



(180) 



The usefulness of this ratio lies in its interesting range of variation: < v < 1/4, with v = 1/4 in 
the case of equal masses, and v — > in the "test-mass" limit for one of the bodies. 

The 3PN and even 3.5PN center-of-mass equations of motion are obtained by replacing in the 
general- frame 3.5PN equations of motion I|l(j8fl the positions and velocities by their center-of-mass 
expressions, applying as usual the order-reduction of all accelerations where necessary. We write 
the relative acceleration in the center-of-mass frame in the form 

f = -^ [il+ A)n< + B ^} + 0^), (181) 
and find 43 that the coefficients A and B are 
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Up to the 2.5PN order the result agrees with the calculation of |152) . The 3.5PN term is issued from 
Refs. [T331 C3H HSU HZU [H2I dU At the 3PN order we have some gauge-dependent logarithms 
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containing a constant r' which is the "logarithmic barycenter" of the two constants r[ and r' 2 : 

lnr' Q = X 1 \nr[+X 2 \nr' 2 . (184) 

The logarithms in Equations i|182l 11831) . together with the constant r' therein, can be removed 
by applying the gauge transformation (|169|) , while still staying within the class of harmonic coor- 
dinates. The resulting modification of the equations of motion will affect only the coefficients of 
the 3PN order in Equations 1)1821 I183|l . let us denote them by ^3pn and /E^pn- The new values 
of these coefficients, say A' 3PN and B' 3PN , obtained after removal of the logarithms by the latter 
harmonic gauge transformation, are then 
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These gauge-transformed coefficients are useful because they do not yield the usual complica- 
tions associated with logarithms. However, they must be handled with care in applications such 
as |158) . since one must ensure that all other quantities in the problem (energy, angular mo- 
mentum, gravitational-wave fluxes, etc.) are defined in the same specific harmonic gauge avoiding 
logarithms. In the following we shall no longer use the coordinate system leading to Equations 1)1851 
I186fl . Therefore all expressions we shall derive below, notably all those concerning the radiation 
field, are valid in the "standard" harmonic coordinate system in which the equations of motion are 
given by Equation fTHH)) or l|IS2| HS3 ) . 



9.4 Equations of motion and energy for circular orbits 

Most inspiralling compact binaries will have been circularized by the time they become visible by 
the detectors LIGO and VIRGO. In the case of orbits that are circular - apart from the gradual 
2.5PN radiation- reaction inspiral - the complicated equations of motion simplify drastically, since 
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we have r — (nv) — 0(l/c 5 ), and the remainder can always be neglected at the 3PN level. In the 
case of circular orbits, up to the 2.5PN order, the relation between center-of-mass variables and 
the relative ones reads |16) 33 
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To display conveniently the successive post-Newtonian corrections, we employ the post-Newtonian 
parameter 

7^ = o(~Y (188) 
rc z V c / 

Notice that there are no corrections of order 1PN in Equations <|187[1 f° r circular orbits; the 
dominant term is of order 2PN, i.e. proportional to 7 2 = C(l/c 4 ). 

The relative acceleration a 1 = a\ — a\ of two bodies moving on a circular orbit at the 3PN 
order is then given by 

32G 3 m^^ o fl_ 
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where x % = y\ — y\ is the relative separation (in harmonic coordinates) and u> denotes the angular 
frequency of the circular motion. The second term in Equation l|189f) . opposite to the velocity 



v l = v\ — v l 2 , is the 2.5PN radiation reaction force (we neglect here its 3.5PN extension), which 
comes from the reduction of the coefficient of 1 /c 5 in Equations (|182l I183|) . The main content of 
the 3PN equations i|189fl is the relation between the frequency u> and the orbital separation r, that 
we find to be given by the generalized version of Kepler's third law |37l 138) : 
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The length scale r' is given in terms of the two gauge-constants r[ and r' 2 by Equation (|184|) . As for 
the energy, it is immediately obtained from the circular-orbit reduction of the general result (|170|) . 
We have 
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(191) 



This expression is that of a physical observable E; however, it depends on the choice of a coordi- 
nate system, as it involves the post-Newtonian parameter 7 defined from the harmonic-coordinate 

33 Actually, in the present computation we do not need the radiation-reaction 2.5PN term in these relations; we 
give it only for completeness. 
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separation ri2. But the numerical value of E should not depend on the choice of a coordinate 
system, so E must admit a frame-invariant expression, the same in all coordinate systems. To find 
it we re-express E with the help of a frequency-related parameter x instead of the post-Newtonian 
parameter 7. Posing 
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we readily obtain from Equation 119Ufl the expression of 7 in terms of x at 3PN order, 
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that we substitute back into Equation (|191|) . making all appropriate post-Newtonian re-expansions. 
As a result, we gladly discover that the logarithms together with their associated gauge constant 



r' have cancelled out. Therefore, our result is 
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For circular orbits one can check that there are no terms of order x 7 ^ 2 in Equation l|194fl . so our 
result for E is actually valid up to the 3.5PN order. 



9.5 The innermost circular orbit (ICO) 

Having in hand the circular-orbit energy, we define the innermost circular orbit (ICO) as the 
minimum, when it exists, of the energy function E(x). Notice that we do not define the ICO as 
a point of dynamical general-relativistic unstability. Hence, we prefer to call this point the ICO 
rather than, strictly speaking, an innermost stable circular orbit or ISCO. A study of the dynamical 
stability of circular binary orbits in the post-Newtonian approximation of general relativity can be 
found in Ref. g3]. 

The previous definition of the ICO is motivated by our comparison with the results of numerical 
relativity. Indeed we shall confront the prediction of the standard (Taylor-based) post-Newtonian 
approach with a recent result of numerical relativity by Gourgoulhon, Grandclement, and Bonaz- 
zola |12(JI 1123] , These authors computed numerically the energy of binary black holes under the 
assumptions of conformal flatness for the spatial metric and of exactly circular orbits. The latter 
restriction is implemented by requiring the existence of an "helical" Killing vector, which is time- 
like inside the light cylinder associated with the circular motion, and space-like outside. In the 
numerical approach [1201 1123] there are no gravitational waves, the field is periodic in time, and 
the gravitational potentials tend to zero at spatial infinity within a restricted model equivalent 
to solving five out of the ten Einstein field equations (the so-called Isenberg- Wilson-Mathews ap- 
proximation; see Ref. |114) for a discussion). Considering an evolutionary sequence of equilibrium 
configurations Refs. [1201 1123] obtained numerically the circular-orbit energy E{u>) and looked for 
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the ICO of binary black holes (see also Refs. [521 11211 |l51j for related calculations of binary neutron 
and strange quark stars). 

Since the numerical calculation [121)1 I123| has been performed in the case of corotating black 
holes, which are spinning with the orbital angular velocity to, we must for the comparison include 
within our post-Newtonian formalism the effects of spins appropriate to two Kerr black holes 
rotating at the orbital rate. The total relativistic mass of the Kerr black hole is given by 34 

M 2 = M? T + -^g-, (195) 

irr 

where S is the spin, related to the usual Kerr parameter by S = Ma, and M- m is the irreducible 
mass given by M- m = \[Aj (47r) (A is the hole's surface area). The angular velocity of the corotating 
black hole is u> — dM/dS hence, from Equation l|195fl . 

c 

(196) 
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Physically this angular velocity is the one of the outgoing photons that remain for ever at the 
location of the light-like horizon. Combining Equations i|195lll9b1l we obtain M and S as functions 
of M; rr and w, 



M = 



2 w 2 ' 



(197) 

s _ 4M 3 r u, 

v /l-4A/ 2 r o; 2 ' 

This is the right thing to do since u> is the basic variable describing each equilibrium configuration 
calculated numerically, and because the irreducible masses are the ones which are held constant 
along the numerical evolutionary sequences in Refs. [12011123] . In the limit of slow rotation we get 

S = Iuj + 0(oj 3 ), (198) 

where / = 4M^ r is the moment of inertia of the black hole. Next the total mass-energy is 

M = M iTI + ^Iuj 2 + (w 4 ) , (199) 

which involves, as we see, the usual kinetic energy of the spin. 

To take into account the spin effects our first task is to replace all the masses entering the 
energy function 1[194[) by their equivalent expressions in terms of u> and the two irreducible masses. 
It is clear that the leading contribution is that of the spin kinetic energy given by Equation (|199[l , 
and it comes from the replacement of the rest mass-energy mc 2 (where m — M\ + M2). From 
Equation ijl99[l this effect is of order us 2 in the case of corotating binaries, which means by compar- 
ison with Equation (|194[) that it is equivalent to an "orbital" effect at the 2PN order (i.e. oc x 2 ). 
Higher-order corrections in Equation i|199[l . which behave at least like w 4 , will correspond to the 
orbital 5PN order at least and are negligible for the present purpose. In addition there will be 
a subdominant contribution, of the order of lu 8 / 3 equivalent to 3PN order, which comes from the 
replacement of the masses into the "Newtonian" part, proportional to a; oc cj 2 / 3 , of the energy E 
(see Equation (|194[) ~). With the 3PN accuracy we do not need to replace the masses that enter 
into the post-Newtonian corrections in E, so in these terms the masses can be considered to be 
the irreducible ones. 



In this section we pose G = 1 = c, and the two individual black hole masses are denoted Mi and M2. 



73 



Our second task is to include the specific relativistic effects due to the spins, namely the spin- 
orbit (SO) interaction and the spin-spin (SS) one. In the case of spins Si and S2 aligned parallel 
to the orbital angular momentum (and right-handed with respect to the sense of motion) the SO 
energy reads 

'4 \ Si /4 Mf \ S 2 ' 



E SO = ~H (mcj) 5 / 3 



3 m 2 / Mf \3m 2 J Mf 



(200) 



Here we are employing the formula given by Kidder et al. [1431 1141) (based on seminal works of 
Barker and O'Connell 0|S]) wn0 have computed the SO contribution and expressed it by means 
of the orbital frequency u>. The derivation of Equation (|200l) in Ref. |143l 1141] takes into account 
the fact that the relation between the orbital separation r (in the harmonic coordinate system) 
and the frequency uj depends on the spins. We immediately infer from Equation l|200|l that in the 
case of corotating black holes the SO effect is equivalent to a 3PN orbital effect and thus must 
be retained with the present accuracy (with this approximation, the masses in Equation (|200|l are 
the irreducible ones). As for the SS interaction (still in the case of spins aligned with the orbital 
angular momentum) it is given by 

E ss =,Hm^^ l . (201) 

The SS effect can be neglected here because it is of order 5PN for corotating systems. Summing up 
all the spin contributions we find that the suplementary energy due to the corotating spins is |23| 

AE COIOt = mc 2 x{{2-Qv)x 2 + (-61s + Wis 2 ) x 3 + 0{x 4 )} , (202) 

where x = (iw) 2 ' 3 . The complete 3PN energy of the corotating binary is finally given by the 
sum of Equations \ 194(1 and l(202|l , in which we must now understand all the masses as being the 
irreducible ones (we no longer indicate the superscript "irr"), which for the comparison with the 
numerical calculation must be assumed to stay constant when the binary evolves. 

The Figure H (issued from Ref. [22]) presents our results for £ico in the case of irrotational 
and corotational binaries. Since AE COIOt , given by Equation l(202|l . is at least of order 2PN, the 
result for lPN corot is the same as for 1PN in the irrotational case; then, obviously, 2PN corot 
takes into account only the leading 2PN corotation effect (i.e. the spin kinetic energy given by 
Equation (I199I) - ). while 3pj\f corot involves also, in particular, the corotational SO coupling at the 
3PN order. In addition we present in Figure^the numerical point obtained by numerical relativity 
under the assumptions of conformal flatness and of helical symmetry |120II123] . As we can see the 
3PN points, and even the 2PN ones, are rather close to the numerical value. The fact that the 
2PN and 3PN values are so close to each other is a good sign of the convergence of the expansion; 
we shall further comment this point in Section IHlfl In fact one might say that the role of the 3PN 
approximation is merely to "confirm" the value already given by the 2PN one (but of course, had 
we not computed the 3PN term, we would not be able to trust very much the 2PN value). As 
expected, the best agreement we obtain is for the 3PN approximation and in the case of corotation, 
i.e. the point 3PN col ° . However, the 1PN approximation is clearly not precise enough, but this is 
not surprising in the highly relativistic regime of the ICO. 

In conclusion, we find that the location of the ICO as computed by numerical relativity, under 
the helical-symmetry and conformal-flatness approximations, is in good agreement with the post- 
Newtonian prediction. See also Ref. [SS] for the results calculated within the effective-one-body 
approach method 16 1) at the 3PN order, which are close to the ones reported in Figure ^ This 
agreement constitutes an appreciable improvement of the previous situation, because the earlier 
estimates of the ICO in post- Newtonian theory |142| and numerical relativity [1771 1§] strongly dis- 
agreed with each other, and do not match with the present 3PN results. The numerical calculation 
of quasi-equilibrium configurations has been since then redone and refined by a number of groups, 
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Figure 1: Results for the binding energy -Eico versus loico i n the equal-mass case (v — l/4j. The 
asterisk marks the result calculated by numerical relativity. The points indicated by 1PN, 2PN ; and 
3PN are computed from the minimum of Equation j!94\ ), and correspond to irrotational binaries. 
The points denoted by ipN corot ; 2PN corot , and 3PN corot come from the minimum of the sum of 
Equations \194\ ) and and describe corotational binaries. 
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for both corotational and irrotational binaries (see in particular Ref. |74] ) . These works confirm 
the previous findings. 

9.6 Accuracy of the post-Newtonian approximation 

In this section we want to assess the validity of the post-Newtonian approximation, and, more 
precisely, to address, and to some extent to answer, the following questions: How accurate is the 
post-Newtonian expansion for describing the dynamics of binary black hole systems? Is the ICO 
of binary black holes, defined by the minimum of the energy function E(u>), accurately determined 
at the highest currently known post-Newtonian order? The latter question is pertinent because 
the ICO represents a point in the late stage of evolution of the binary which is very relativistic 
(orbital velocities of the order of 50% of the speed of light). How well does the 3PN approximation 
as compared with the prediction provided by numerical relativity (see Section I9.5J1 ? What is 
the validity of the various post-Newtonian resummation techniques |92l 1931 1601 l6*T] which aim at 
"boosting" the convergence of the standard post-Newtonian approximation? 

The previous questions are interesting but difficult to settle down rigorously. Indeed the very 
essence of an approximation is to cope with our ignorance of the higher-order terms in some 
expansion, but the higher-order terms are precisely the ones which would be needed for a satisfying 
answer to these problems. So we shall be able to give only some educated guesses and/or plausible 
answers, that we cannot justify rigorously, but which seem very likely from the standard point 
of view on the post-Newtonian theory, in particular that the successive orders of approximation 
get smaller and smaller as they should (in average), with maybe only few accidents occuring 
at high orders where a particular approximation would be abnormally large with respect to the 
lower-order ones. Admittedly, in addition, our faith in the estimation we shall give regarding the 
accuracy of the 3PN order for instance, comes from the historical perspective, thanks to the many 
successes achieved in the past by the post-Newtonian approximation when confronting the theory 
and observations. It is indeed beyond question, from our past experience, that the post-Newtonian 
method does work. 

Establishing the post-Newtonian expansion rigorously has been the subject of numerous math- 
ematical oriented works, see, e.g., [1841 11851 1186) . In the present section we shall simply look 
(much more modestly) at what can be said by inspection of the explicit post-Newtonian coeffi- 
cients which have been computed so far. Basically, the point we would like to emphasize 35 is that 
the post-Newtonian approximation, in standard form (without using the resummation techniques 
advocated in Refs. 1601 is able to located the ICO of two black holes, in the case of 

comparable masses (mi ~ m^), with a very good accuracy. At first sight this statement is rather 
surprising, because the dynamics of two black holes at the point of the ICO is so relativistic. In- 
deed one sometimes hears about the "bad convergence" , or the "fundamental breakdown" , of the 
post-Newtonian series in the regime of the ICO. However our estimates do show that the 3PN 
approximation is good in this regime, for comparable masses, and we have already confirmed this 
by the remarkable agreement with the numerical calculations, as detailed in Section 19.51 

Let us center our discussion on the post-Newtonian expression of the circular-orbit energy (1 194(1 , 
developed to the 3PN order, which is of the form 

2 

E(x) = -Hl£ {1 + ai (v) x + a 2 { V ) x 2 + a 3 (u) x 3 + 0(x 4 )} . (203) 
The first term, proportional to a;, is the Newtonian term, and then we have many post-Newtonian 

35 We are following the discussion in Ref. 1241 . Note that the arguments of this section are rather biased toward 
the author's own work 1231 1241 . 



7G 



corrections, the coefficients of which are known up to 3PN order [Hfil IT37I 1351 1371 ETFI 11113] : 

3 2/ 
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For the discussion it is helpful to keep the Hadamard rcgularization ambiguity parameter A present 
in the 3PN coefficient a^{v). Recall from Section l^2*l that this parameter was introduced in Refs. [23 
EH] and is equivalent to the parameter astatic of Refs. |136lll37j . We already gave in Equation ljl33J) 
the relation linking them, 

^- 3 1987 fon<\ 

A - --(Static ~ 3080' ( 205 ) 

Before its actual computation in general relativity, it has been argued in Ref. that the 
numerical value of astatic could be ~ —9, because for such a value some different resummation 
techniques, when they are implemented at the 3PN order, give approximately the same result for 
the ICO. Even more, it was suggested [51] that w s t a tic might be precisely equal to w* tatic , with 

47 41 

Static = - y + = -9.34 • • • . (206) 

However, as reviewed in Sections 18 . 21 and 18 .31 the computations performed using dimensional regu- 
larization, within the ADM-Hamiltonian formalism |96) and harmonic-coordinate approaches |30) . 
and the independent computation of Refs. [1301 1129"] . have settled the value of this parameter in 
general relativity to be 

1987 

Static =0 A = -^- ( 20? ) 

We note that this result is quite different from w* tatic , Equation (|206|) . This already suggests that 
different resummation techniques, namely Pade approximants |92l 1931 194| and effective-one-body 
methods |6UI 1611 l9"l] . which are designed to "accelerate" the convergence of the post-Newtonian 
series, do not in fact converge toward the same exact solution (or, at least, not as fast as expected). 

In the limiting case v — ► 0, the expression (|203l 1204^1 reduces to the 3PN approximation of the 
energy for a test particle in the Schwarzschild background, 



E Sch (x) =^c 2 



1 - 2x 



vT - 3.x 



(208) 



The minimum of that function or Schwarzschild ICO occurs at xfco = 1/6, and we have Ef^Q = 
lie 2 (a/8/9 — l). We know that the Schwarzschild ICO is also an innermost stable circular orbit 
or ISCO, i.e. it corresponds to a point of dynamical unstability. Another important feature of 
Equation (|208|l is the singularity at the value ^fightring = 1/3 which corresponds to the famous 
circular orbit of photons in the Schwarzschild metric ( "light-ring" singularity) . This orbit can also 
be viewed as the last unstable circular orbit. We can check that the post-Newtonian coefficients 



a 



Sch 



a„(0) corresponding to Equation H208|) are given by 



3 sch_ 3"(2n-l)!!(2n-l) 
" 2»(n+l)! ■ 1 J 

They increase with n by roughly a factor 3 at each order. This is simply the consequence of the fact 
that the radius of convergence of the post-Newtonian series is given by the Schwarzschild light-ring 
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singularity at the value 1/3. We may therefore recover the light-ring orbit by investigating the 
limit 

a Sc \ 1 

~ZSoh = Q = flight ring- (210) 



Let us now discuss a few order-of-magnitude estimates. At the location of the ICO we have 
found (see Figure^in Section l9~5jl that the frequency-related parameter x defined by Equation l|192[l 
is approximately of the order of x ~ (0.1) 2 / 3 ~ 20% for equal masses. Therefore, we might a priori 
expect that the contribution of the 1PN approximation to the energy at the ICO should be of that 
order. For the present discussion we take the pessimistic view that the order of magnitude of an 
approximation represents also the order of magnitude of the higher-order terms which are neglected. 
We see that the 1PN approximation should yield a rather poor estimate of the "exact" result, but 
this is quite normal at this very relativistic point where the orbital velocity is v/c ~ x 1 / 2 ~ 50%. 
By the same argument we infer that the 2PN approximation should do much better, with fractional 
errors of the order of x 2 ~ 5%, while 3PN will be even better, with the accuracy a; 3 ~ 1%. 

Now the previous estimate makes sense only if the numerical values of the post-Newtonian 
coefficients in Equations l|204|l stay roughly of the order of one. If this is not the case, and 
if the coefficients increase dangerously with the post-Newtonian order n, one sees that the post- 
Newtonian approximation might in fact be very bad. It has often been emphasized in the litterature 
(see, e.g., Refs. |Z3 EHSI E2J ) that in the test -mass limit v — > the post-Newtonian series converges 
slowly, so the post-Newtonian approximation is not very good in the regime of the ICO. Indeed 
we have seen that when v = the radius of convergence of the series is 1/3 (not so far from 
x ico = 1/6); an d that accordingly the post-Newtonian coefficients increase by a factor ~ 3 at 
each order. So it is perfectly correct to say that in the case of test particles in the Schwarzschild 
background the post-Newtonian approximation is to be carried out to a high order in order to 
locate the turning point of the ICO. 

What happens when the two masses are comparable (u = 4)? It is clear that the accuracy of the 
post-Newtonian approximation depends crucially on how rapidly the post-Newtonian coefficients 
increase with n. We have seen that in the case of the Schwarzschild metric the latter increase is 
in turn related to the existence of a light-ring orbit. For continuing the discussion we shall say 
that the relativistic interaction between two bodies of comparable masses is "Schwarzschild-like" 
if the post-Newtonian coefficients a n (j) increase when n — > +oo. If this is the case this signals the 
existence of something like a light-ring singularity which could be interpreted as the deformation, 
when the mass ratio v is "turned on", of the Schwarzschild light-ring orbit. By analogy with 
Equation l|210|l we can estimate the location of this "pseudo-light-ring" orbit by 

" 1,1 flight ring (y) with n = 3. (211) 



a n (v) 

Here n = 3 is the highest known post-Newtonian order. If the two-body problem is "Schwarzschild- 
like" then the right-hand side of Equation (|211|) is small (say around 1/3), the post-Newtonian 
coefficients typically increase with n, and most likely it should be difficult to get a reliable estimate 
by post-Newtonian methods of the location of the ICO. So we ask: Is the gravitational interaction 
between two comparable masses Schwarzschild-like? 

In Table ^ we present the values of the coefficients a n (v) in the test-mass limit v = (see 
Equation l|209|l for their analytic expression), and in the equal-mass case v = j when the ambiguity 
parameter takes the "uncorrect" value w* tatic defined by Equation (|206|> . and the correct one 
^static = predicted by general relativity. When is = we clearly see the expected increase of the 
coefficients by roughly a factor 3 at each step. Now, when v = \ and ^static = ^static we notice 
that the coefficients increase approximately in the same manner as in the test- mass case v = 0. 
This indicates that the gravitational interaction in the case of w* tatic looks like that in a one-body 
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Table 1 : Numerical values of the sequence of coefficients of the post- Newtonian series composing 
the energy function E(x) as given by Equations ^2031 \204j )- 



problem. The associated pseudo-light-ring singularity is estimated using Equation i|211[l as 

flight ring Static J ~ °- 32 - ( 212 ) 

The pseudo-light-ring orbit seems to be a very small deformation of the Schwarzschild light-ring 
orbit given by Equation l|210|l . In this Schwarzschild-like situation, we should not expect the 
post-Newtonian series to be very accurate. 

Now in the case v = \ but when the ambiguity parameter takes the correct value w sta tic = 0, we 
see that the 3PN coefficient 013(7;) is of the order of —1 instead of being ~ —10. This suggests, unless 
3PN happens to be quite accidental, that the post-Newtonian coefficients in general relativity do 
not increase very much with n. This is an interesting finding because it indicates that the actual 
two-body interaction in general relativity is not Schwarzschild-like. There does not seem to exist 
something like a light-ring orbit which would be a deformation of the Schwarzschild one. Applying 
Equation l|211|l we obtain as an estimate of the "light ring" , 

^lightring f-,GR^ ~ 2.88. (213) 

It is clear that if we believe the correctness of this estimate we must conclude that there is in 
fact no notion of a light-ring orbit in the real two-body problem. Or, one might say (pictorially 
speaking) that the light-ring orbit gets hidden inside the horizon of the final black hole formed by 
coalescence. Furthermore, if we apply Equation (|211() using the 2PN approximation n — 2 instead 
of the 3PN one n = 3, we get the value ~ 0.28 instead of Equation (|213|) . So at the 2PN order the 
metric seems to admit a light ring, while at the 3PN order it apparently does not admit any. This 
erratic behaviour reinforces our idea that it is meaningless (with our present 3PN-based knowledge, 
and until fuller information is available) to assume the existence of a light-ring singularity when 
the masses are equal. 

It is impossible of course to be thoroughly confident about the validity of the previous state- 
ment because we know only the coefficients up to 3PN order. Any tentative conclusion based on 
3PN can be "falsified" when we obtain the next 4PN order. Nevertheless, we feel that the mere 
fact that aj,{j) — —0.97 in Tabled is sufficient to motivate our conclusion that the gravitational 
field generated by two bodies is more complicated than the Schwarzschild space-time. This ap- 
praisal should look cogent to relativists and is in accordance with the author's respectfulness of 
the complexity of the Einstein field equations. 

We want next to comment on a possible implication of our conclusion as regards the so-called 
post-Newtonian resummation techniques, i.e. Pade approximants 92 , 93 , 94 , which aim at "boost- 
ing" the convergence of the post-Newtonian series in the pre-coalescence stage, and the effective- 
one-body (EOB) method [HOII^IM]! which attempts at describing the late stage of the coalescence 
of two black holes. These techniques are based on the idea that the gravitational two-body inter- 
action is a "deformation" - with v < j being the deformation parameter - of the Schwarzschild 
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space-time. The Pade approximants are valuable tools for giving accurate representations of func- 
tions having some singularities. In the problem at hands they would be justified if the "exact" 
expression of the energy (whose 3PN expansion is given by Equations (|203l I204[) ) would admit 
a singularity at some reasonable value of x (e.g., < 0.5). In the Schwarzschild case, for which 
Equation H210[l holds, the Pade series converges rapidly The Pade constructed only from the 
2PN approximation of the energy - keeping only af ch and of 1311 - already coincide with the exact 
result given by Equation H208J1 . On the other hand, the EOB method maps the post-Newtonian 
two-body dynamics (at the 2PN or 3PN orders) on the geodesic motion on some effective metric 
which happens to be a i/-deformation of the Schwarzschild space-time. In the EOB method the 
effective metric looks like Schwarzschild by definition, and we might of course expect the two-body 
interaction to own the main Schwarzschild-like features. 

Our comment is that the validity of these post-Newtonian resummation techniques does not 
seem to be compatible with the value w s t a ti c = 0, which suggests that the two-body interaction 
in general relativity is not Schwarzschild-like. This doubt is confirmed by the finding of Ref. |94| 
(already alluded to above) that in the case of the wrong ambiguity parameter w* tatic ~ —9.34 
the Pade approximants and the EOB method at the 3PN order give the same result for the ICO. 
From the previous discussion we see that this agreement is to be expected because a deformed 
light-ring singularity seems to exist with that value ^g tatic . By contrast, in the case of general 
relativity, (jJ sta tic = 0, the Pade and EOB methods give quite different results (cf. the Figure 2 
in |94|). Another confirmation comes from the light-ring singularity which is determined from the 
Pade approximants at the 2PN order (see Equation (3.22) in [H2]) as 



This value is rather close to Equation l|212|l but strongly disagrees with Equation l|213|l . Our 
explanation is that the Pade series converges toward a theory having astatic — ^l tat - ic ; such a 
theory is different from general relativity. 

Finally we come to the good news that, if really the post-Newtonian coefficients when v — \ stay 
of the order of one (or minus one) as it seems to, this means that the standard post-Newtonian 
approach, based on the standard Taylor approximants, is probably very accurate. The post- 
Newtonian series is likely to "converge well", with a "convergence radius" of the order of one 36 . 
Hence the order-of-magnitude estimate we proposed at the beginning of this section is probably 
correct. In particular the 3PN order should be close to the "exact" solution for comparable masses 
even in the regime of the ICO. 



36 Actually, the post-Newtonian series could be only asymptotic (hence divergent), but nevertheless it should give 
excellent results provided that the series is truncated near some optimal order of approximation. In this discussion 
we assume that the 3PN order is not too far from that optimum. 




(214) 
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10 Gravitational Waves from Compact Binaries 



We pointed out that the 3PN equations of motion, Equations (|189l I190)) . are merely Newtonian 
as regards the radiative aspects of the problem, because with that precision the radiation reaction 
force is at the lowest 2.5PN order. A solution would be to extend the precision of the equations of 
motion so as to include the full relative 3PN or 3.5PN precision into the radiation reaction force, 
but, needless to say, the equations of motion up to the 5.5PN or 6PN order are quite impossible 
to derive with the present technology. The much better alternative solution is to apply the wave- 
generation formalism described in Part [fl and to determine by its means the work done by the 
radiation reaction force directly as a total energy flux at future null infinity. In this approach, we 
replace the knowledge of the higher-order radiation reaction force by the computation of the total 
flux C, and we apply the energy balance equation as in the test of the P of the binary pulsar (see 
Equations |@1EJ|): 

f~* (215, 

Therefore, the result (|194l) that we found for the 3.5PN binary's center-of-mass energy E constitutes 
only "half" of the solution of the problem. The second "half" consists of finding the rate of decrease 
dE/dt, which by the balance equation is nothing but finding the total gravitational- wave flux C at 
the 3.5PN order. Because the orbit of inspiralling binaries is circular, the balance equation for the 
energy is sufficient (no need of a balance equation for the angular momentum). This all sounds 
perfect, but it is important to realize that we shall use Equation (|215|) at the very high 3.5PN 
order, at which order there are no proofs (following from first principles in general relativity) that 
the equation is correct, despite its physically obvious character. Nevertheless, Equation i|215[l has 
been checked to be valid, both in the cases of point-particle binaries [13311134] and extended weakly 
self-gravitating fluids ^JE|, a t the 1PN order and even at 1.5PN (the 1.5PN approximation is 
especially important for this check because it contains the first wave tails). 

Obtaining C can be divided into two equally important steps: (1) the computation of the source 
multipole moments II and Jl of the compact binary and (2) the control and determination of the 
tails and related non-linear effects occuring in the relation between the binary's source moments 
and the radiative ones Ul and (cf. the general formalism of Part [J. 

10.1 The binary's multipole moments 

The general expressions of the source multipole moments given by Theorem Equations l|85[l . are 
first to be worked out explicitly for general fluid systems at the 3PN order. For this computation 
one uses the formula (|91|) . and we insert the 3PN metric coefficients (in harmonic coordinates) 
expressed in Equations (|115l) by means of the retarded- type elementary potentials l|117l 11181 1119|) . 
Then we specialize each of the (quite numerous) terms to the case of point-particle binaries by 
inserting, for the matter stress-energy tensor T a @ , the standard expression made out of Dirac 
delta- functions. The infinite self-field of point-particles is removed by means of the Hadamard 
regularization; and dimensional regularization is used to compute the few ambiguity parameters 
(see Section [S}. This computation has been performed in at the 1PN order, and in at the 
2PN order; we report below the most accurate 3PN results obtained in Refs. [451 1441 15T1 132) . 

The difficult part of the analysis is to find the closed-form expressions, fully explicit in terms of 
the particle's positions and velocities, of many non-linear integrals. We refer to @H] for full details; 
nevertheless, let us give a few examples of the type of technical formulas that are employed in this 
calculation. Typically we have to compute some integrals like 




(216) 
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where n = |x — yi| and = |x — y2 1 . When n > —3 and p > —3, this integral is perfectly 
well-defined (recall that the finite part TV deals with the bound at infinity). When n < —3 or 
p < —3, our basic ansatz is that we apply the definition of the Hadamard partie finie provided by 
Equation l|124[l . Two examples of closed- form formulas that we get, which do not necessitate the 
Hadamard partie finie, are (quadrupole case I = 2) 
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We denote for example y[ 1 ^ = yfy 3 i (and ri2 = r|yi — y2 1); the constant tq is the one pertaining 



to the finite-part process (see Equation H36[T 
location of the particle 1 is 
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One example where the integral diverges at the 
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where s\ is the Hadamard-regularization constant introduced in Equation (|124|) 37 . 

The crucial input of the computation of the flux at the 3PN order is the mass quadrupole 
moment Iy , since this moment necessitates the full 3PN precision. The result of Ref . @H1 for this 
moment (in the case of circular orbits) is 
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where we pose Xi — x l = y\ 2 and — v l = v\ 2 - The third term is the 2.5PN radiation-reaction 
term, which does not contribute to the energy flux for circular orbits. The two important coeffi- 
cients are A and B, whose expressions through 3PN order are 
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These expressions are valid in harmonic coordinates via the post-Newtonian parameter 7 given by 
Equation (|188|) . As we see, there are two types of logarithms in the moment: One type involves the 
length scale r' related by Equation (|184|) to the two gauge constants r[ and r' 2 present in the 3PN 

37 Whcn computing the gravitational-wave flux in Ref. 1451 we preferred to call the Hadamard-regularization 
constants ui and 112, in order to distinguish them from the constants s± and S2 that were used in our previous 
computation of the equations of motion in Ref. 1381 . Indeed these regularization constants need not neccessarily be 
the same when employed in different contexts. 
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equations of motion; the other type contains the different length scale ro coming from the general 
formalism of Part[I]- indeed, recall that there is a J-V operator in front of the source multipole 
moments in Theorem |SJ As we know, that r is pure gauge; it will disappear from our physical 
results at the end. On the other hand, we have remarked that the multipole expansion outside 
a general post-Newtonian source is actually free of ro, since the ro's present in the two terms of 
Equation l|67|l cancel out. We shall indeed find that the constants ro present in Equations H220[l 
are compensated by similar constants coming from the non-linear wave "tails of tails" . Finally, 
the constants £ , k, and £ are the Hadamard-regularization ambiguity parameters which take the 
values (|136() . 

Besides the 3PN mass quadrupole (l2~T31 . we need also the mass octupole moment Iy^ and 
current quadrupole moment , both of them at the 2PN order; these are given by |45| 
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Also needed are the 1PN mass 2 4 -pole, 1PN current 2 3 -pole (octupole), Newtonian mass 2 5 -pole 
and Newtonian current 2 4 -pole: 
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These results permit one to control what can be called the "instantaneous" part, say Anst, °f 
the total energy flux, by which we mean that part of the flux that is generated solely by the source 
multipole moments, i.e. not counting the "non-instantaneous" tail integrals. The instantaneous 
flux is defined by the replacement into the general expression of C given by Equation l|60|l of all the 
radiative moments \Jl and Vl by the corresponding (Zth time derivatives of the) source moments 
II and J^. Actually, we prefer to define Anst by means of the intermediate moments Mi and Sl- 
Up to the 3.5PN order we have 
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The time derivatives of the source moments (|219l 12201 12211 1222|) are computed by means of the 
circular-orbit equations of motion given by Equation l|189l I190|) : then we substitute them into 
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Equation (J223I> 3S - The net result is 
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The Newtonian approximation, £n = (32c 5 /5G) v 2 ^j b , is the prediction of the Einstein quadrupole 
formula Q, as computed by Landau and Lifchitz |150) . In Equation l|224[) . we have replaced 
the Hadamard regularization ambiguity parameters A and 8 arising at the 3PN order by their 
values and (fH7) . 



10.2 Contribution of wave tails 

To the "instantaneous" part of the flux, we must add the contribution of non-linear multipole 
interactions contained in the relationship between the source and radiative moments. The needed 
material has already been provided in Equations i|§71 19"5|) . Up to the 3.5PN level we have the 
dominant quadratic-order tails, the cubic-order tails or tails of tails, and the non-linear memory 
integral. We shall see that the tails play a crucial role in the predicted signal of compact bina- 
ries. By contrast, the non-linear memory effect, given by the integral inside the 2.5PN term in 
Equation l|97|l . does not contribute to the gravitational- wave energy flux before the 4PN order in 
the case of circular-orbit binaries (essentially because the memory integral is actually "instanta- 
neous" in the flux), and therefore has rather poor observational consequences for future detections 
of inspiralling compact binaries. We split the energy flux into the different terms 
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where Anst has just been found in Equation (|224|) ; £taii is made of the quadratic (multipolar) 
tail integrals in Equation l|98|l : £(taii) 2 is the square of the quadrupole tail in Equation (|97|) : and 
£taii(taii) is the quadrupole tail of tail in Equation l|97|) . We find that £taii contributes at the 
half-integer 1.5PN, 2.5PN, and 3.5PN orders, while both £( ta ii) 2 and £taii(taii) appear only at the 
3PN order. It is quite remarkable that so small an effect as a "tail of tail" should be relevant to 
the present computation, which is aimed at preparing the ground for forthcoming experiments. 

The results follow from the reduction to the case of circular compact binaries of the general 
formulas i|97l I98fl . in which we make use of the explicit expressions for the source moments of 
compact binaries as found in Section ll0.ll Without going into accessory details (see Ref. P3), let 
us point out that following the general formalism of Part [H the total mass M in front of the tail 
integrals is the ADM mass of the binary which is given by the sum of the rest masses, m = mi + m2 
(which is the one appearing in the 7-parameter, Equation l|188[) ). plus some relativistic corrections. 
At the 2PN relative order needed here to compute the tail integrals we have 



M 



1 ~ h + i (7 ~ v) 72 



■04 



(226) 



'For circular orbits there is no difference at this order between 1^, J^, and M^, S^. 
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Let us give the two basic technical formulas needed when carrying out this reduction: 



+00 
dr In: 



+00 
dr In 2 t e~ 



= (C + lnff), 

a 



In a) 



(227) 



where a € C and C = 0.577 • • • denotes the Euler constant |122) . The tail integrals are evaluated 
thanks to these formulas for a fixed (non-decaying) circular orbit. Indeed it can be shown jSO] that 
the "remote-past" contribution to the tail integrals is negligible; the errors due to the fact that the 
orbit actually spirals in by gravitational radiation do not affect the signal before the 4PN order. We 
then find, for the quadratic tail term stricto sensu, the 1.5PN, 2.5PN, and 3.5PN contributions 39 
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For the sum of squared tails and cubic tails of tails at 3PN, we get 
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By comparing Equations (|224l) and (|229|) we observe that the constants ro cleanly cancel out. 
Adding together all these contributions we obtain 
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The gauge constant Tq has not yet disappeared because the post-Newtonian expansion is still 
parametrized by 7 instead of the frequency-related parameter x defined by Equation Ijl92|l ~ just 
as for E when it was given by Equation Ijl91|l . After substituting the expression 7(2;) given by 
Equation (|193(l . we find that r' does cancel as well. Because the relation 7(2;) is issued from the 
equations of motion, the latter cancellation represents an interesting test of the consistency of the 
two computations, in harmonic coordinates, of the 3PN multipole moments and the 3PN equations 
of motion. At long last we obtain our end result: 
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39 All formulas incorporate the changes in some equations following the published Errata (2005) to the works 1161 
rT9ll45ll40l l4). 
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In the test-mass limit v — ► for one of the bodies, we recover exactly the result following from 
linear black- hole perturbations obtained by Tagoshi and Sasaki |202| . In particular, the rational 
fraction 6643739519/69854400 comes out exactly the same as in black-hole perturbations. On the 
other hand, the ambiguity parameters A and 9 are part of the rational fraction —134543/7776, 
belonging to the coefficient of the term at 3PN order proportional to v (hence this coefficient 
cannot be computed by linear black hole perturbations) 40 . 



10.3 Orbital phase evolution 

We shall now deduce the laws of variation with time of the orbital frequency and phase of an 
inspiralling compact binary from the energy balance equation 1)215(1 . The center-of-mass energy E 
is given by Equation 119411 and the total flux C by Equation (|231l) . For convenience we adopt the 
dimcnsionlcss time variable 41 
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where t c denotes the instant of coalescence, at which the frequency tends to infinity (evidently, the 
post-Newtonian method breaks down well before this point). We transform the balance equation 
into an ordinary differential equation for the parameter x, which is immediately integrated with 
the result 
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The orbital phase is defined as the angle 4>, oriented in the sense of the motion, between the 
separation of the two bodies and the direction of the ascending node Af within the plane of the 
sky, namely the point on the orbit at which the bodies cross the plane of the sky moving toward 

40 Generalizing the flux formula 12311 to point masses moving on quasi elliptic orbits dates back to the work of 
Peters and Mathews 11751 at Newtonian order. The result was obtained in 12141 14lj| at 1PN order, and then further 
extended by Gopakumar and Iyer 11191 up to 2PN order using an explicit quasi-Keplerian representation of the 
motion 99 194 . No complete result at 3PN order is yet available. 

41 Notice the "strange" post-Newtonian order of this time variable: = C(c+ 8 ). 
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the detector. We have dcj>/dt — to, which translates, with our notation, into d(f>/dO = —5/u 
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from which we determine 
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where Qq is a constant of integration that can be fixed by the initial conditions when the wave 
frequency enters the detector's bandwidth. Finally we want also to dispose of the important 
expression of the phase in terms of the frequency x. For this we get 
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where xq is another constant of integration. With the formula (|235(l the orbital phase is complete 
up to the 3.5PN order. The effects due to the spins of the particles, i.e. the spin-orbit (SO) 
coupling arising at the 1.5PN order for maximally rotating compact bodies and the spin-spin (SS) 
coupling at the 2PN order, can be added if necessary; they are known up to the 2.5PN order 
included [T4lfl H2T1 IPd51 I2UT1 ITTU1 125] . On the other hand, the contribution of the quadrupole 
moments of the compact objects, which are induced by tidal effects, is expected to come only at 
the 5PN order (see Equation ©). 

As a rough estimate of the relative importance of the various post-Newtonian terms, let us give 
in Table their contributions to the accumulated number of gravitational-wave cycles J\f in the 
bandwidth of the LIGO and VIRGO detectors (see also Table I in Ref. 35 for the contributions 
of the SO and SS effects). Note that such an estimate is only indicative, because a full treatment 
would require the knowledge of the detector's power spectral density of noise, and a complete 
simulation of the parameter estimation using matched filtering |79l 11811 1149) . We define AT by 
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The frequency of the signal at the entrance of the bandwidth is the seismic cut-off frequency 
/seismic of ground-based detectors; the terminal frequency /isco is assumed for simplicity's sake 
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to be given by the Schwarzschild innermost stable circular orbit. Here / = lo/tt = 2/P is the 
signal frequency at the dominant harmonics (twice the orbital frequency). As we see in Tabled 
with the 3PN or 3.5PN approximations we reach an acceptable level of, say, a few cycles, that 
roughly corresponds to the demand which was made by data-analysists in the case of neutron-star 
binaries |77l 1781 1751 ASHI 1591 15%] . Indeed, the above estimation suggests that the neglected 4PN 
terms will yield some systematic errors that are, at most, of the same order of magnitude, i.e. a 
few cycles, and perhaps much less (see also the discussion in Section l9~Hf . 
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1.5PN (dominant tail) 
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2PN 
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3PN 
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2.2 


3.5PN 
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Table 2: Contributions of post-Newtonian orders to the accumulated number of gravitational-wave 
cycles M (defined by Equation $236)) ) in the bandwidth of VIRGO and LIGO detectors. Neutron 
stars have mass 1.4 M©, and black holes 10 M©. The entry frequency is / se ismic = 10 Hz, and the 
terminal frequency is /isco = c 3 /(6 3 / 2 7rGm). 



10.4 The two polarization waveforms 

The theoretical templates of the compact binary inspiral follow from insertion of the previous 
solutions for the 3.5PN-accurate orbital frequency and phase into the binary's two polarization 
waveforms h+ and h x ■ We shall include in /i+ and h x all the harmonics, besides the dominant one 
at twice the orbital frequency, up to the 2.5PN order, as they have been calculated in Refs. [4l>ll4"|. 
The polarization waveforms are defined with respect to two polarization vectors p = (pi) and 

q = fe), 
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where p and q are chosen to lie along the major and minor axis, respectively, of the projection 
onto the plane of the sky of the circular orbit, with p oriented toward the ascending node Af. To 
the 2PN order we have 

= {<l + * 1/2 <x 2) + + > + + + o (1) } . 

(238) 

The post- Newtonian terms are ordered by means of the frequency- related variable x. They depend 
on the binary's 3.5PN-accurate phase <f> through the auxiliary phase variable 
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where M = m [l — v^j2 + (l/c 4 )] is the ADM mass (cf. Equation 1|226|) ). and where ojq is a con- 
stant frequency that can conveniently be chosen to be the entry frequency of a laser-interferometric 
detector (say wo/ir = 10 Hz). For the plus polarization we have 42 
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For the cross polarization, we have 
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42 We neglect the non-linear memory (DC) term present in the Newtonian plus polarization See Wiseman 

and Will 12191 and Arun et al. [1] for the computation of this term. 
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3 Sin 

~-^ s i c i {I 3 + 10 ln 2 ] cos -0 + 5?r sin ^ - 9 [7 - 10 ln(3/2)l cos 30; - 45?r sin 30} . 

20 m 

We use the shorthands Cj = cos i and Sj = sin i for the cosine and sine of the inclination angle i 
between the direction of the detector as seen from the binary's center-of-mass, and the normal to 
the orbital plane (we always suppose that the normal is right-handed with respect to the sense of 
motion, so that < i < ir). Finally, the more recent calculation of the 2.5PN order in Ref. 0] is 
reported here: 



H 



(5/2) 



5m 

Si COS 01 



1771 1667 



5120 5120 



217 
9216' 



-7T COS 2lp 



19 

y 



/681 13 




35 


^256 + 768 




768 


( 3451 


673 




V 9216 


3072 


-**-§«*- 


-( 


16 

~y 



4 + 



9216 
1 



2304 



1 



9216 
14 



-Si — cos dip 
m 



3537 


22977 


1024 


5120 
' 23829 




1280 


+ ^ 2 


/29127 


[ 5120 



cf- 



15309 
5120 £ 



3072 
729 



5120 



cos 40 



- — (1 + C?)^(1 



5m , 
+Si — cos 50 
m 



108125 
9216 4 
/8125 
V 256 



5529 , 
1280 Cq 
27267 
5120 C 

- dv) 

40625 



7749 4 
1280 C ' 1 ' 
1647 4 
" 5120 Ci 



729 , 
1280 C ' 

2187 
" 5120 C 



9216 1 
40625 
" 2304 ' 



83125 4 
9216 Ct 
48125 



15625 6 
9216 q 
15625 



2304 



2304 



119375 40625 



-d 



■ cos 70 



117649 



46080 



+ sin 20 
+s?(l + c 2 )sin40 



9216 3072 
4(1 + ^X1-4^ + 3^) 
96 8 , 



44375 
9216 ' 



15625 
3072 ' 



14 2 7 4 

y^+r^ 



28 



56 

y 



32 In 2 



5 1 
1193 

~yr 



5 

-32 In 2 



(240) 



H 



(5/2) 



6 2 



+d cos 20 



+CjSj cos 40; 

Sm . 

+SiCi sm0 

m 



22 



154 94 



112 64, „ 
H n 2 

5 3 



5 

1193 
15 



64 ln2 



913 



1891 



7680 11520 



4608 



90 



1165 235 



384 576 



1152 



v 2 \- 



1301 301 



+7TCi sin 2-0 



34 

t 



3 Ci 



4608 2304 



1536 



-SiCi sin dtp 

m 



12501 12069 



2560 1280 



1701 
2560 ' 



-v 



19581 7821 



-v 



640 320 ' 
2 / 18903 11403 
V 2560 ~ 1280 ' 



1701 , 
640 C '' 
5103 

" 2560' 



-sfci sin 4-0 
5m 

s^i sin Dip 



32tt 



(1 - 31/) 



101875 6875 



4608 256 1 
( 66875 44375 



21875 
4608 C 



V 1152 



576 



21875 
1152 1 



100625 83125 



4608 



2304 



21875 
1536 ' 



s^Ci smtip 

171 



117649 
23040 



(1 - 4i/ + 3z/ 2 



(241) 



The practical implementation of the theoretical templates in the data analysis of detectors 
follows the standard matched filtering technique. The raw output of the detector o(t) consists 
of the superposition of the real gravitational wave signal h Iea x(t) and of noise n(t). The noise is 
assumed to be a stationary Gaussian random variable, with zero expectation value, and with (sup- 
posedly known) frequency-dependent power spectral density S n (u). The experimenters construct 
the correlation between o(t) and a filter q(t), i.e. 



dt' o(t')q(t + t'), 



(242) 



and divide c(t) by the square root of its variance, or correlation noise. The expectation value of 
this ratio defines the filtered signal-to- noise ratio (SNR). Looking for the useful signal h rea \(t) in 
the detector's output o(t), the experimenters adopt for the filter 



q(uj) 
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(243) 



where q(oj) and h(u) are the Fourier transforms of q(t) and of the theoretically computed template 
h(t). By the matched filtering theorem, the filter (|243|) maximizes the SNR if hit) — h Tea }(t). The 
maximum SNR is then the best achievable with a linear filter. In practice, because of systematic 
errors in the theoretical modelling, the template hit) will not exactly match the real signal h Ioa \(t), 
but if the template is to constitute a realistic representation of nature the errors will be small. 
This is of course the motivation for computing high order post-Newtonian templates, in order to 
reduce as much as possible the systematic errors due to the unknown post-Newtonian remainder. 

To conclude, the use of theoretical templates based on the preceding 2.5PN wave forms, and 
having their frequency evolution built in via the 3.5PN phase evolution (|234l l235Tl . should yield 
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some accurate detection and measurement of the binary signals. Interestingly, it should also permit 
some new tests of general relativity, because we have the possibility of checking that the observed 
signals do obey each of the terms of the phasing formulas l|234l I235|l . e.g., those associated with 
the specific non-linear tails, exactly as they are predicted by Einstein's theory 07|0;8 5 . Indeed, 
we don't know of any other physical systems for which it would be possible to perform such tests. 
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